Walmart Weekly Sales Forecasting: Model Comparison

Author

Machine Learning UTEC Homework

1 Objective

This report compares six forecasting approaches for weekly Walmart sales:

  1. Local Linear Anomaly Model (Phase 1).
  2. Bayesian Structural AR (Phase 2, full data).
  3. Elastic Net anomaly baseline (Phase 3).
  4. Random Forest anomaly baseline (Phase 4).
  5. ETS anomaly baseline (Phase 5).
  6. AdaBoost anomaly baseline (Phase 7).

Additionally, we include a Phase 1-lite diagnostic variant in the final comparison to show the impact of removing the extra exogenous variables from Local Linear.

Evaluation uses forward-chaining validation and weighted MAE (WMAE).

2 Notebook Training Cells (Scratch-First, Function-Based)

These cells are self-contained and intended for notebook execution from scratch.

Key design choices: - No argparse and no sys.argv runner. - Hyperparameters are defined in notebook cells. - Workflow per phase is: run baseline config -> tune on CV -> save best config -> run final model. - outputs/*/config*.json are written as outputs, not read as inputs.

2.1 Shared Setup

from pathlib import Path
import json

import numpy as np
import pandas as pd

PROJECT_ROOT = Path(".").resolve()


def _coerce_markdown_series(series: pd.Series) -> pd.Series:
    if series.dtype == bool:
        return series.astype(float)
    as_str = series.astype(str).str.strip().str.upper()
    mapped = as_str.map({"TRUE": 1.0, "FALSE": 0.0, "NA": np.nan, "NAN": np.nan, "": np.nan})
    numeric = pd.to_numeric(series, errors="coerce")
    if numeric.notna().sum() == 0 and mapped.notna().sum() > 0:
        return mapped.astype(float)
    out = numeric.astype(float)
    mask = out.isna() & mapped.notna()
    out.loc[mask] = mapped.loc[mask]
    return out


def ensure_prepared_features(root: Path = PROJECT_ROOT) -> tuple[Path, Path]:
    train_feat_path = root / "train_feat.parquet"
    test_feat_path = root / "test_feat.parquet"
    if train_feat_path.exists() and test_feat_path.exists():
        return train_feat_path, test_feat_path

    train_csv = root / "train.csv"
    test_csv = root / "test.csv"
    features_csv = root / "features.csv"
    stores_csv = root / "stores.csv"
    required = [train_csv, test_csv, features_csv, stores_csv]
    missing = [str(p) for p in required if not p.exists()]
    if missing:
        raise FileNotFoundError(f"Missing required raw files: {missing}")

    train = pd.read_csv(train_csv)
    test = pd.read_csv(test_csv)
    features = pd.read_csv(features_csv)
    stores = pd.read_csv(stores_csv)

    for frame in (train, test, features):
        frame["Date"] = pd.to_datetime(frame["Date"], errors="coerce")

    for frame in (train, test, features):
        if "IsHoliday" in frame.columns:
            frame["IsHoliday"] = (
                frame["IsHoliday"]
                .astype(str)
                .str.upper()
                .map({"TRUE": True, "FALSE": False})
                .fillna(frame["IsHoliday"])
            )

    for col in ["MarkDown1", "MarkDown2", "MarkDown3", "MarkDown4", "MarkDown5"]:
        if col in features.columns:
            features[col] = _coerce_markdown_series(features[col])

    train_feat = train.merge(features, on=["Store", "Date", "IsHoliday"], how="left")
    train_feat = train_feat.merge(stores, on=["Store"], how="left")
    test_feat = test.merge(features, on=["Store", "Date", "IsHoliday"], how="left")
    test_feat = test_feat.merge(stores, on=["Store"], how="left")

    train_feat.to_parquet(train_feat_path, index=False)
    test_feat.to_parquet(test_feat_path, index=False)
    return train_feat_path, test_feat_path


def mean_row_from_metrics(metrics_df: pd.DataFrame) -> dict:
    if metrics_df.empty:
        raise ValueError("metrics_df is empty")
    mask = metrics_df["fold"].astype(str) == "mean"
    if mask.any():
        return metrics_df.loc[mask].iloc[0].to_dict()
    return metrics_df.iloc[-1].to_dict()


train_feat_path, test_feat_path = ensure_prepared_features(PROJECT_ROOT)
print(f"Using train features: {train_feat_path}")
print(f"Using test features:  {test_feat_path}")

2.2 Phase 1 Core Definitions

import json
from collections import deque
from dataclasses import dataclass
from pathlib import Path

import numpy as np
import pandas as pd

TARGET_COL = "sales_anom"
LAG_ORDERS = (1, 2)
LAG_FEATURE_COLS = [f"sales_anom_lag{lag}" for lag in LAG_ORDERS]
FULL_EXOG_FEATURE_COLS = [
    "temp_anom",
    "fuel_anom",
    "MarkDown1",
    "MarkDown2",
    "MarkDown3",
    "MarkDown4",
    "MarkDown5",
    "CPI",
    "Unemployment",
]
LITE_EXOG_FEATURE_COLS = [
    "temp_anom",
    "fuel_anom",
]
CORE_FEATURE_COLS = [
    "sales_anom_lag1",
    "sales_anom_lag2",
    "is_holiday_int",
]
INTERACTION_FEATURE_COLS = [
    "holiday_x_lag1",
]
WEEK_PERIOD = 53

@dataclass(slots=True)
class LocalLinearConfig:
    kernel: str
    bandwidth: int
    min_samples: int
    ridge: float
    coef_clip: float
    anom_clip_scale: float
    use_interactions: bool
    feature_mode: str

@dataclass(slots=True)
class FittedLocalLinear:
    local_params: dict[tuple[str, int], np.ndarray]
    available_weeks: dict[str, np.ndarray]
    series_fallback: dict[str, np.ndarray]
    global_fallback: np.ndarray
    series_bounds: dict[str, tuple[float, float]]
    global_bounds: tuple[float, float]
    feature_cols: list[str]
    feature_means: np.ndarray
    feature_stds: np.ndarray

def get_feature_cols(config: LocalLinearConfig) -> list[str]:
    cols = list(get_exog_feature_cols(config)) + list(CORE_FEATURE_COLS)
    if config.use_interactions:
        cols.extend(INTERACTION_FEATURE_COLS)
    return cols

def get_exog_feature_cols(config: LocalLinearConfig) -> list[str]:
    if config.feature_mode == "lite":
        return list(LITE_EXOG_FEATURE_COLS)
    return list(FULL_EXOG_FEATURE_COLS)

def weighted_ridge_fit(X: np.ndarray, y: np.ndarray, weights: np.ndarray, ridge: float) -> np.ndarray:
    if X.size == 0:
        return np.zeros(X.shape[1] + 1, dtype=float)

    design = np.column_stack([np.ones(X.shape[0], dtype=float), X])
    sqrt_w = np.sqrt(np.clip(weights, 0.0, None))
    Xw = design * sqrt_w[:, None]
    yw = y * sqrt_w

    reg = np.eye(Xw.shape[1], dtype=float) * ridge
    reg[0, 0] = 0.0

    xtx = Xw.T @ Xw + reg
    xty = Xw.T @ yw
    beta = np.linalg.pinv(xtx) @ xty
    return beta.astype(float, copy=False)

def regularize_beta(beta: np.ndarray, coef_clip: float) -> np.ndarray:
    out = beta.copy()
    out[1:] = np.clip(out[1:], -coef_clip, coef_clip)
    return out

def cyclic_week_distance(weeks: np.ndarray, target_week: int) -> np.ndarray:
    delta = np.abs(weeks - target_week)
    return np.minimum(delta, WEEK_PERIOD - delta)

def kernel_weights(distance: np.ndarray, bandwidth: int, kernel: str) -> np.ndarray:
    if kernel == "tricube":
        scaled = np.clip(distance / float(bandwidth), 0.0, 1.0)
        weights = (1.0 - scaled**3) ** 3
        weights[distance > bandwidth] = 0.0
        return weights
    if kernel == "exponential":
        return np.exp(-0.5 * (distance / float(bandwidth)) ** 2)
    raise ValueError(f"Unsupported kernel: {kernel}")

def robust_bounds(values: np.ndarray, clip_scale: float) -> tuple[float, float]:
    if values.size == 0:
        return (-1.0, 1.0)
    q01, q99 = np.quantile(values, [0.01, 0.99])
    std = float(np.std(values))
    span = max(float(q99 - q01), std, 1.0)
    lower = float(q01 - clip_scale * span)
    upper = float(q99 + clip_scale * span)
    if lower >= upper:
        center = float(np.median(values))
        lower, upper = center - 1.0, center + 1.0
    return (lower, upper)

def fit_standard_scaler(X: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    means = np.mean(X, axis=0)
    stds = np.std(X, axis=0)
    bad = (~np.isfinite(stds)) | (stds < 1e-8)
    stds[bad] = 1.0
    means[~np.isfinite(means)] = 0.0
    return means.astype(float, copy=False), stds.astype(float, copy=False)

def apply_standard_scaler(X: np.ndarray, means: np.ndarray, stds: np.ndarray) -> np.ndarray:
    return (X - means) / stds

def load_data(train_path: Path, test_path: Path, max_series: int | None) -> tuple[pd.DataFrame, pd.DataFrame]:
    train = pd.read_parquet(train_path)
    test = pd.read_parquet(test_path)

    for frame in (train, test):
        frame["Date"] = pd.to_datetime(frame["Date"], errors="coerce")
        frame["Store"] = frame["Store"].astype(int)
        frame["Dept"] = frame["Dept"].astype(int)
        frame["series_id"] = frame["Store"].astype(str) + "_" + frame["Dept"].astype(str)
        frame["week_of_year"] = frame["Date"].dt.isocalendar().week.astype(int)
        frame["is_holiday_int"] = frame["IsHoliday"].astype(int)

    if max_series is not None:
        keep_ids = (
            train["series_id"]
            .drop_duplicates()
            .sort_values()
            .head(max_series)
            .tolist()
        )
        train = train[train["series_id"].isin(keep_ids)].copy()
        test = test[test["series_id"].isin(keep_ids)].copy()

    train = train.sort_values(["Date", "Store", "Dept"]).reset_index(drop=True)
    test = test.sort_values(["Date", "Store", "Dept"]).reset_index(drop=True)
    return train, test

def build_forward_folds(
    train_df: pd.DataFrame, n_folds: int, val_weeks: int
) -> list[tuple[int, np.ndarray, np.ndarray]]:
    unique_dates = np.array(sorted(train_df["Date"].dropna().unique()))
    initial_train = len(unique_dates) - (n_folds * val_weeks)
    if initial_train <= 0:
        raise ValueError(
            f"Not enough unique weeks ({len(unique_dates)}) for n_folds={n_folds}, val_weeks={val_weeks}."
        )

    folds: list[tuple[int, np.ndarray, np.ndarray]] = []
    for fold in range(n_folds):
        train_end = initial_train + fold * val_weeks
        val_end = train_end + val_weeks
        train_dates = unique_dates[:train_end]
        val_dates = unique_dates[train_end:val_end]
        folds.append((fold + 1, train_dates, val_dates))
    return folds

def compute_climatology(train_df: pd.DataFrame) -> dict[str, pd.DataFrame | dict[str, float]]:
    series = (
        train_df.groupby(["series_id", "week_of_year"], as_index=False)[
            ["Weekly_Sales", "Temperature", "Fuel_Price"]
        ]
        .median()
        .rename(
            columns={
                "Weekly_Sales": "clim_sales_series",
                "Temperature": "clim_temp_series",
                "Fuel_Price": "clim_fuel_series",
            }
        )
    )
    store = (
        train_df.groupby(["Store", "week_of_year"], as_index=False)[
            ["Weekly_Sales", "Temperature", "Fuel_Price"]
        ]
        .median()
        .rename(
            columns={
                "Weekly_Sales": "clim_sales_store",
                "Temperature": "clim_temp_store",
                "Fuel_Price": "clim_fuel_store",
            }
        )
    )
    global_week = (
        train_df.groupby("week_of_year", as_index=False)[["Weekly_Sales", "Temperature", "Fuel_Price"]]
        .median()
        .rename(
            columns={
                "Weekly_Sales": "clim_sales_global",
                "Temperature": "clim_temp_global",
                "Fuel_Price": "clim_fuel_global",
            }
        )
    )
    overall = {
        "sales": float(train_df["Weekly_Sales"].median()),
        "temp": float(train_df["Temperature"].median()),
        "fuel": float(train_df["Fuel_Price"].median()),
    }
    return {"series": series, "store": store, "global_week": global_week, "overall": overall}

def attach_climatology(
    df: pd.DataFrame, climatology: dict[str, pd.DataFrame | dict[str, float]]
) -> pd.DataFrame:
    out = df.merge(
        climatology["series"], on=["series_id", "week_of_year"], how="left"  # type: ignore[arg-type]
    )
    out = out.merge(climatology["store"], on=["Store", "week_of_year"], how="left")  # type: ignore[arg-type]
    out = out.merge(climatology["global_week"], on=["week_of_year"], how="left")  # type: ignore[arg-type]

    overall = climatology["overall"]  # type: ignore[assignment]
    out["clim_sales"] = (
        out["clim_sales_series"].fillna(out["clim_sales_store"]).fillna(out["clim_sales_global"]).fillna(overall["sales"])
    )
    out["clim_temp"] = (
        out["clim_temp_series"].fillna(out["clim_temp_store"]).fillna(out["clim_temp_global"]).fillna(overall["temp"])
    )
    out["clim_fuel"] = (
        out["clim_fuel_series"].fillna(out["clim_fuel_store"]).fillna(out["clim_fuel_global"]).fillna(overall["fuel"])
    )

    out["temp_anom"] = out["Temperature"] - out["clim_temp"]
    out["fuel_anom"] = out["Fuel_Price"] - out["clim_fuel"]
    if "Weekly_Sales" in out.columns:
        out["sales_anom"] = out["Weekly_Sales"] - out["clim_sales"]

    drop_cols = [
        "clim_sales_series",
        "clim_temp_series",
        "clim_fuel_series",
        "clim_sales_store",
        "clim_temp_store",
        "clim_fuel_store",
        "clim_sales_global",
        "clim_temp_global",
        "clim_fuel_global",
    ]
    return out.drop(columns=drop_cols, errors="ignore")

def add_anomaly_lags(df: pd.DataFrame) -> pd.DataFrame:
    out = df.sort_values(["series_id", "Date"]).copy()
    grp = out.groupby("series_id")["sales_anom"]
    for lag in LAG_ORDERS:
        out[f"sales_anom_lag{lag}"] = grp.shift(lag)
    out["holiday_x_lag1"] = out["is_holiday_int"] * out["sales_anom_lag1"]
    return out

def fit_local_models(train_df: pd.DataFrame, config: LocalLinearConfig) -> FittedLocalLinear:
    feature_cols = get_feature_cols(config)
    valid = train_df.copy()
    # Keep early rows by treating unavailable long-lag anomalies as neutral (0).
    for c in LAG_FEATURE_COLS + INTERACTION_FEATURE_COLS:
        if c in valid.columns:
            valid[c] = valid[c].fillna(0.0)
    for c in get_exog_feature_cols(config):
        if c in valid.columns and valid[c].isna().any():
            median = float(valid[c].median())
            if not np.isfinite(median):
                median = 0.0
            valid[c] = valid[c].fillna(median)
    valid = valid.dropna(subset=feature_cols + [TARGET_COL]).copy()
    if valid.empty:
        zero_beta = np.zeros(len(feature_cols) + 1, dtype=float)
        return FittedLocalLinear(
            {},
            {},
            {},
            zero_beta,
            {},
            (-1.0, 1.0),
            feature_cols,
            np.zeros(len(feature_cols), dtype=float),
            np.ones(len(feature_cols), dtype=float),
        )

    all_X_raw = valid[feature_cols].to_numpy(dtype=float)
    scaler_means, scaler_stds = fit_standard_scaler(all_X_raw)
    all_X = apply_standard_scaler(all_X_raw, scaler_means, scaler_stds)
    all_y = valid[TARGET_COL].to_numpy(dtype=float)
    global_fallback = regularize_beta(
        weighted_ridge_fit(
        all_X, all_y, np.ones(all_y.shape[0], dtype=float), ridge=config.ridge
        ),
        coef_clip=config.coef_clip,
    )

    anom_values = train_df[TARGET_COL].dropna().to_numpy(dtype=float)
    global_bounds = robust_bounds(anom_values, clip_scale=config.anom_clip_scale)
    series_bounds: dict[str, tuple[float, float]] = {}

    local_params: dict[tuple[str, int], np.ndarray] = {}
    available_weeks: dict[str, np.ndarray] = {}
    series_fallback: dict[str, np.ndarray] = {}

    grouped = valid.groupby("series_id", sort=False)
    for i, (series_id, group) in enumerate(grouped, start=1):
        g = group.sort_values("Date")
        X_raw = g[feature_cols].to_numpy(dtype=float)
        X = apply_standard_scaler(X_raw, scaler_means, scaler_stds)
        y = g[TARGET_COL].to_numpy(dtype=float)
        weeks = g["week_of_year"].to_numpy(dtype=int)
        series_bounds[series_id] = robust_bounds(y, clip_scale=config.anom_clip_scale)

        if y.shape[0] >= 3:
            series_fallback[series_id] = regularize_beta(
                weighted_ridge_fit(X, y, np.ones(y.shape[0], dtype=float), ridge=config.ridge),
                coef_clip=config.coef_clip,
            )
        else:
            series_fallback[series_id] = global_fallback

        learned_weeks: list[int] = []
        for target_week in range(1, WEEK_PERIOD + 1):
            distance = cyclic_week_distance(weeks, target_week)
            weights = kernel_weights(distance, config.bandwidth, config.kernel)
            active = weights > 0.0
            if int(active.sum()) < config.min_samples:
                continue
            beta = regularize_beta(
                weighted_ridge_fit(X[active], y[active], weights[active], ridge=config.ridge),
                coef_clip=config.coef_clip,
            )
            local_params[(series_id, target_week)] = beta
            learned_weeks.append(target_week)

        if learned_weeks:
            available_weeks[series_id] = np.array(learned_weeks, dtype=int)

        if i % 300 == 0:
            print(f"  fitted models for {i} series...")

    return FittedLocalLinear(
        local_params,
        available_weeks,
        series_fallback,
        global_fallback,
        series_bounds,
        global_bounds,
        feature_cols,
        scaler_means,
        scaler_stds,
    )

def _lag_from_history(hist: deque[float], lag: int) -> float:
    if len(hist) >= lag:
        return float(hist[-lag])
    return 0.0

def initialize_history(train_df: pd.DataFrame) -> dict[str, deque[float]]:
    history: dict[str, deque[float]] = {}
    max_lag = max(LAG_ORDERS)
    ordered = train_df.sort_values(["series_id", "Date"])
    for series_id, group in ordered.groupby("series_id", sort=False):
        values = group["sales_anom"].dropna().to_numpy(dtype=float)
        history[series_id] = deque(values.tolist(), maxlen=max_lag)
    return history

def pick_beta(model: FittedLocalLinear, series_id: str, week_of_year: int) -> np.ndarray:
    direct = model.local_params.get((series_id, week_of_year))
    if direct is not None:
        return direct

    learned = model.available_weeks.get(series_id)
    if learned is not None and learned.size > 0:
        delta = np.abs(learned - week_of_year)
        cyclic = np.minimum(delta, WEEK_PERIOD - delta)
        nearest_week = int(learned[np.argmin(cyclic)])
        return model.local_params[(series_id, nearest_week)]

    series_fb = model.series_fallback.get(series_id)
    if series_fb is not None:
        return series_fb
    return model.global_fallback

def _build_feature_vector(row, hist: deque[float], feature_cols: list[str]) -> np.ndarray:
    def safe_float(value: float) -> float:
        value = float(value)
        if np.isfinite(value):
            return value
        return 0.0

    lag_vals = {
        "sales_anom_lag1": _lag_from_history(hist, 1),
        "sales_anom_lag2": _lag_from_history(hist, 2),
    }
    feat = {
        "temp_anom": safe_float(row.temp_anom),
        "fuel_anom": safe_float(row.fuel_anom),
        "MarkDown1": safe_float(row.MarkDown1),
        "MarkDown2": safe_float(row.MarkDown2),
        "MarkDown3": safe_float(row.MarkDown3),
        "MarkDown4": safe_float(row.MarkDown4),
        "MarkDown5": safe_float(row.MarkDown5),
        "CPI": safe_float(row.CPI),
        "Unemployment": safe_float(row.Unemployment),
        "is_holiday_int": safe_float(row.is_holiday_int),
        **lag_vals,
        "holiday_x_lag1": safe_float(row.is_holiday_int) * lag_vals["sales_anom_lag1"],
    }
    return np.array([feat[c] for c in feature_cols], dtype=float)

def recursive_predict(
    future_df: pd.DataFrame, model: FittedLocalLinear, history: dict[str, deque[float]]
) -> pd.DataFrame:
    ordered = future_df.sort_values(["Date", "series_id"]).copy()
    preds: list[tuple[int, float, float]] = []
    max_lag = max(LAG_ORDERS)

    for row in ordered.itertuples(index=False):
        series_id = row.series_id
        hist = history.get(series_id)
        if hist is None:
            hist = deque(maxlen=max_lag)
            history[series_id] = hist

        beta = pick_beta(model, series_id, int(row.week_of_year))
        x_raw = _build_feature_vector(row, hist, model.feature_cols)
        x = apply_standard_scaler(x_raw, model.feature_means, model.feature_stds)
        pred_anom = float(beta[0] + np.dot(beta[1:], x))
        lower, upper = model.series_bounds.get(series_id, model.global_bounds)
        pred_anom = float(np.clip(pred_anom, lower, upper))
        pred_sales = float(pred_anom + row.clim_sales)

        preds.append((int(row.row_id), pred_sales, pred_anom))
        hist.append(pred_anom)

    pred_df = pd.DataFrame(preds, columns=["row_id", "pred_sales", "pred_anom"])
    return future_df.merge(pred_df, on="row_id", how="left")

def wmae(y_true: np.ndarray, y_pred: np.ndarray, holiday: np.ndarray) -> float:
    weights = np.where(holiday, 5.0, 1.0)
    return float(np.average(np.abs(y_true - y_pred), weights=weights))

def summarize_metrics(df: pd.DataFrame) -> dict[str, float]:
    y_true = df["Weekly_Sales"].to_numpy(dtype=float)
    y_pred = df["pred_sales"].to_numpy(dtype=float)
    holiday = df["IsHoliday"].to_numpy(dtype=bool)

    metrics = {
        "wmae": wmae(y_true, y_pred, holiday),
        "mae": float(np.mean(np.abs(y_true - y_pred))),
        "rmse": float(np.sqrt(np.mean((y_true - y_pred) ** 2))),
        "n_rows": float(df.shape[0]),
    }
    if "sales_anom" in df.columns and "pred_anom" in df.columns:
        metrics["anom_mae"] = float(np.mean(np.abs(df["sales_anom"] - df["pred_anom"])))
    return metrics

def run_cv(
    train_df: pd.DataFrame, config: LocalLinearConfig, n_folds: int, val_weeks: int
) -> tuple[pd.DataFrame, pd.DataFrame]:
    folds = build_forward_folds(train_df, n_folds=n_folds, val_weeks=val_weeks)
    metrics_rows: list[dict[str, float | int | str]] = []
    oof_parts: list[pd.DataFrame] = []

    for fold_id, train_dates, val_dates in folds:
        print(f"\nFold {fold_id}: train={len(train_dates)} weeks, val={len(val_dates)} weeks")

        fold_train = train_df[train_df["Date"].isin(train_dates)].copy()
        fold_val = train_df[train_df["Date"].isin(val_dates)].copy()
        fold_val["row_id"] = np.arange(fold_val.shape[0], dtype=int)

        climatology = compute_climatology(fold_train)
        fold_train = add_anomaly_lags(attach_climatology(fold_train, climatology))

        model = fit_local_models(fold_train, config)
        history = initialize_history(fold_train)

        fold_val = attach_climatology(fold_val, climatology)
        fold_val_pred = recursive_predict(fold_val, model, history)

        fold_metrics = summarize_metrics(fold_val_pred)
        fold_metrics["fold"] = fold_id
        fold_metrics["train_start"] = str(pd.Timestamp(train_dates[0]).date())
        fold_metrics["train_end"] = str(pd.Timestamp(train_dates[-1]).date())
        fold_metrics["val_start"] = str(pd.Timestamp(val_dates[0]).date())
        fold_metrics["val_end"] = str(pd.Timestamp(val_dates[-1]).date())
        metrics_rows.append(fold_metrics)

        print(
            f"  wmae={fold_metrics['wmae']:.4f} mae={fold_metrics['mae']:.4f} "
            f"rmse={fold_metrics['rmse']:.4f}"
        )

        keep_cols = [
            "Store",
            "Dept",
            "Date",
            "IsHoliday",
            "Weekly_Sales",
            "sales_anom",
            "pred_sales",
            "pred_anom",
        ]
        part = fold_val_pred[keep_cols].copy()
        part["fold"] = fold_id
        oof_parts.append(part)

    metrics_df = pd.DataFrame(metrics_rows)
    if not metrics_df.empty:
        mean_row = {"fold": "mean"}
        for col in ["wmae", "mae", "rmse", "anom_mae", "n_rows"]:
            if col in metrics_df.columns:
                mean_row[col] = float(metrics_df[col].mean())
        metrics_df = pd.concat([metrics_df, pd.DataFrame([mean_row])], ignore_index=True)

    oof_df = pd.concat(oof_parts, ignore_index=True) if oof_parts else pd.DataFrame()
    return metrics_df, oof_df

def fit_full_and_predict_test(
    train_df: pd.DataFrame, test_df: pd.DataFrame, config: LocalLinearConfig
) -> pd.DataFrame:
    train_full = train_df.copy()
    test_full = test_df.copy()
    test_full["row_id"] = np.arange(test_full.shape[0], dtype=int)

    climatology = compute_climatology(train_full)
    train_full = add_anomaly_lags(attach_climatology(train_full, climatology))
    model = fit_local_models(train_full, config)
    history = initialize_history(train_full)

    test_full = attach_climatology(test_full, climatology)
    pred = recursive_predict(test_full, model, history)
    return pred

def ensure_dirs(base: Path) -> dict[str, Path]:
    paths = {
        "root": base,
        "metrics": base / "metrics",
        "preds": base / "preds",
        "submissions": base / "submissions",
    }
    for p in paths.values():
        p.mkdir(parents=True, exist_ok=True)
    return paths

def make_submission(pred_df: pd.DataFrame) -> pd.DataFrame:
    out = pred_df.copy()
    out["Id"] = (
        out["Store"].astype(int).astype(str)
        + "_"
        + out["Dept"].astype(int).astype(str)
        + "_"
        + out["Date"].dt.strftime("%Y-%m-%d")
    )
    submission = out[["Id", "pred_sales"]].rename(columns={"pred_sales": "Weekly_Sales"})
    return submission


def phase1_default_config() -> dict:
    return {
        "train_path": "train_feat.parquet",
        "test_path": "test_feat.parquet",
        "output_dir": "outputs",
        "kernel": "tricube",
        "bandwidth": 6,
        "min_samples": 16,
        "ridge": 1e-4,
        "coef_clip": 6.0,
        "anom_clip_scale": 2.0,
        "n_folds": 4,
        "val_weeks": 13,
        "max_series": None,
        "use_interactions": False,
        "feature_mode": "full",
    }


def run_phase1(
    config: dict | None = None,
    _LocalLinearConfig=LocalLinearConfig,
    _ensure_dirs=ensure_dirs,
    _load_data=load_data,
    _run_cv=run_cv,
    _fit_full_and_predict_test=fit_full_and_predict_test,
    _make_submission=make_submission,
    _get_exog_feature_cols=get_exog_feature_cols,
    _get_feature_cols=get_feature_cols,
) -> dict[str, Path]:
    cfg = phase1_default_config()
    if config:
        cfg.update(config)

    ll_cfg = _LocalLinearConfig(
        kernel=cfg["kernel"],
        bandwidth=int(cfg["bandwidth"]),
        min_samples=int(cfg["min_samples"]),
        ridge=float(cfg["ridge"]),
        coef_clip=float(cfg["coef_clip"]),
        anom_clip_scale=float(cfg["anom_clip_scale"]),
        use_interactions=bool(cfg["use_interactions"]),
        feature_mode=str(cfg["feature_mode"]),
    )

    out = _ensure_dirs(Path(cfg["output_dir"]))
    train_df, test_df = _load_data(
        Path(cfg["train_path"]),
        Path(cfg["test_path"]),
        max_series=cfg.get("max_series"),
    )

    metrics_df, oof_df = _run_cv(
        train_df,
        config=ll_cfg,
        n_folds=int(cfg["n_folds"]),
        val_weeks=int(cfg["val_weeks"]),
    )

    metrics_path = out["metrics"] / "local_linear_cv.csv"
    oof_path = out["preds"] / "local_linear_oof.parquet"
    metrics_df.to_csv(metrics_path, index=False)
    oof_df.to_parquet(oof_path, index=False)

    test_pred = _fit_full_and_predict_test(train_df, test_df, config=ll_cfg)
    test_pred_path = out["preds"] / "local_linear_test.parquet"
    test_pred.to_parquet(test_pred_path, index=False)

    submission = _make_submission(test_pred)
    submission_path = out["submissions"] / "local_linear_submission.csv"
    submission.to_csv(submission_path, index=False)

    config_payload = {
        **cfg,
        "lags": list(LAG_ORDERS),
        "interaction_cols": INTERACTION_FEATURE_COLS if cfg["use_interactions"] else [],
        "exogenous_features": _get_exog_feature_cols(ll_cfg),
        "feature_cols": _get_feature_cols(ll_cfg),
        "standard_scaler": True,
    }
    config_path = out["metrics"] / "local_linear_config.json"
    config_path.write_text(json.dumps(config_payload, indent=2), encoding="utf-8")

    return {
        "metrics": metrics_path,
        "oof": oof_path,
        "test": test_pred_path,
        "submission": submission_path,
        "config": config_path,
    }


def tune_phase1(
    base_config: dict,
    search_space: list[dict],
    search_output_dir: str = "outputs/search/phase1",
    _LocalLinearConfig=LocalLinearConfig,
    _load_data=load_data,
    _run_cv=run_cv,
) -> tuple[pd.DataFrame, dict]:
    cfg0 = phase1_default_config()
    cfg0.update(base_config)

    train_df, _ = _load_data(
        Path(cfg0["train_path"]),
        Path(cfg0["test_path"]),
        max_series=cfg0.get("max_series"),
    )

    rows: list[dict] = []
    best_cfg: dict | None = None
    best_wmae = float("inf")

    for idx, cand in enumerate(search_space, start=1):
        cfg = {**cfg0, **cand}
        ll_cfg = _LocalLinearConfig(
            kernel=cfg["kernel"],
            bandwidth=int(cfg["bandwidth"]),
            min_samples=int(cfg["min_samples"]),
            ridge=float(cfg["ridge"]),
            coef_clip=float(cfg["coef_clip"]),
            anom_clip_scale=float(cfg["anom_clip_scale"]),
            use_interactions=bool(cfg["use_interactions"]),
            feature_mode=str(cfg["feature_mode"]),
        )

        metrics_df, _ = _run_cv(
            train_df,
            config=ll_cfg,
            n_folds=int(cfg["n_folds"]),
            val_weeks=int(cfg["val_weeks"]),
        )
        mean_row = mean_row_from_metrics(metrics_df)
        wmae_val = float(mean_row["wmae"])
        rec = {
            "candidate": idx,
            **cand,
            "wmae": wmae_val,
            "mae": float(mean_row.get("mae", np.nan)),
            "rmse": float(mean_row.get("rmse", np.nan)),
        }
        rows.append(rec)

        if wmae_val < best_wmae:
            best_wmae = wmae_val
            best_cfg = cfg

    results_df = pd.DataFrame(rows).sort_values("wmae").reset_index(drop=True)
    out_dir = Path(search_output_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    results_df.to_csv(out_dir / "search_results.csv", index=False)
    if best_cfg is None:
        raise RuntimeError("Phase 1 tuning did not evaluate any candidates")
    (out_dir / "best_config.json").write_text(json.dumps(best_cfg, indent=2), encoding="utf-8")
    return results_df, best_cfg

2.3 Baseline Utility Definitions

import time
from collections import deque
from pathlib import Path
from typing import Callable

import numpy as np
import pandas as pd

TARGET_COL = "sales_anom"
LAG_ORDERS = (1, 2)
EXOG_FEATURE_COLS = [
    "temp_anom",
    "fuel_anom",
    "MarkDown1",
    "MarkDown2",
    "MarkDown3",
    "MarkDown4",
    "MarkDown5",
    "CPI",
    "Unemployment",
]
ANOM_FEATURE_COLS = [
    "lag1",
    "lag2",
    "week_of_year",
    "month",
    "is_holiday_int",
    *EXOG_FEATURE_COLS,
]

def load_data(train_path: Path, test_path: Path, max_series: int | None = None) -> tuple[pd.DataFrame, pd.DataFrame]:
    train = pd.read_parquet(train_path).copy()
    test = pd.read_parquet(test_path).copy()

    for frame in (train, test):
        frame["Date"] = pd.to_datetime(frame["Date"], errors="coerce")
        frame["Store"] = frame["Store"].astype(int)
        frame["Dept"] = frame["Dept"].astype(int)
        frame["series_id"] = frame["Store"].astype(str) + "_" + frame["Dept"].astype(str)
        frame["week_of_year"] = frame["Date"].dt.isocalendar().week.astype(int)
        frame["month"] = frame["Date"].dt.month.astype(int)
        frame["is_holiday_int"] = frame["IsHoliday"].astype(int)

    if max_series is not None:
        keep_ids = train["series_id"].drop_duplicates().sort_values().head(max_series).tolist()
        train = train[train["series_id"].isin(keep_ids)].copy()
        test = test[test["series_id"].isin(keep_ids)].copy()

    train = train.sort_values(["Date", "Store", "Dept"]).reset_index(drop=True)
    test = test.sort_values(["Date", "Store", "Dept"]).reset_index(drop=True)
    return train, test

def build_forward_folds(train_df: pd.DataFrame, n_folds: int, val_weeks: int) -> list[tuple[int, np.ndarray, np.ndarray]]:
    unique_dates = np.array(sorted(train_df["Date"].dropna().unique()))
    initial_train = len(unique_dates) - (n_folds * val_weeks)
    if initial_train <= 0:
        raise ValueError(
            f"Not enough unique weeks ({len(unique_dates)}) for n_folds={n_folds}, val_weeks={val_weeks}."
        )

    folds: list[tuple[int, np.ndarray, np.ndarray]] = []
    for fold in range(n_folds):
        train_end = initial_train + fold * val_weeks
        val_end = train_end + val_weeks
        train_dates = unique_dates[:train_end]
        val_dates = unique_dates[train_end:val_end]
        folds.append((fold + 1, train_dates, val_dates))
    return folds

def ensure_output_dir(path: Path) -> Path:
    path.mkdir(parents=True, exist_ok=True)
    return path

def wmae(y_true: np.ndarray, y_pred: np.ndarray, holiday: np.ndarray) -> float:
    weights = np.where(holiday, 5.0, 1.0)
    return float(np.average(np.abs(y_true - y_pred), weights=weights))

def summarize_metrics(df: pd.DataFrame) -> dict[str, float]:
    y_true = df["Weekly_Sales"].to_numpy(dtype=float)
    y_pred = df["pred_sales"].to_numpy(dtype=float)
    holiday = df["IsHoliday"].to_numpy(dtype=bool)
    metrics = {
        "wmae": wmae(y_true, y_pred, holiday),
        "mae": float(np.mean(np.abs(y_true - y_pred))),
        "rmse": float(np.sqrt(np.mean((y_true - y_pred) ** 2))),
        "n_rows": float(df.shape[0]),
    }
    if "sales_anom" in df.columns and "pred_anom" in df.columns:
        metrics["anom_mae"] = float(np.mean(np.abs(df["sales_anom"] - df["pred_anom"])))
    return metrics

def _get_lag_value(hist: deque[float], lag: int) -> float:
    if len(hist) >= lag:
        return float(hist[-lag])
    return np.nan

def _robust_bounds(values: np.ndarray, clip_scale: float = 2.0) -> tuple[float, float]:
    arr = np.asarray(values, dtype=float)
    arr = arr[np.isfinite(arr)]
    if arr.size == 0:
        return (-1.0, 1.0)
    q01, q99 = np.quantile(arr, [0.01, 0.99])
    std = float(np.std(arr))
    span = max(float(q99 - q01), std, 1.0)
    lower = float(q01 - clip_scale * span)
    upper = float(q99 + clip_scale * span)
    if lower >= upper:
        center = float(np.median(arr))
        lower, upper = center - 1.0, center + 1.0
    return (lower, upper)

def initialize_target_history(train_df: pd.DataFrame, max_lag: int) -> dict[str, deque[float]]:
    history: dict[str, deque[float]] = {}
    ordered = train_df.sort_values(["series_id", "Date"])
    for series_id, group in ordered.groupby("series_id", sort=False):
        values = group[TARGET_COL].dropna().to_numpy(dtype=float)
        history[series_id] = deque(values.tolist(), maxlen=max_lag)
    return history

def make_supervised_lag_frame(
    train_df: pd.DataFrame,
    lag_orders: tuple[int, ...] = LAG_ORDERS,
    feature_cols: list[str] = ANOM_FEATURE_COLS,
) -> tuple[pd.DataFrame, np.ndarray]:
    out = train_df.sort_values(["series_id", "Date"]).copy()
    grouped = out.groupby("series_id")[TARGET_COL]
    for lag in lag_orders:
        out[f"lag{lag}"] = grouped.shift(lag)

    X = out[feature_cols].copy()
    y = out[TARGET_COL].to_numpy(dtype=float)
    valid_mask = np.isfinite(y) & np.all(np.isfinite(X.to_numpy(dtype=float)), axis=1)
    return X.loc[valid_mask].reset_index(drop=True), y[valid_mask]

def recursive_predict_tabular_anomaly(
    future_df: pd.DataFrame,
    model,
    history: dict[str, deque[float]],
    lag_orders: tuple[int, ...] = LAG_ORDERS,
    feature_cols: list[str] = ANOM_FEATURE_COLS,
    lower: float = -np.inf,
    upper: float = np.inf,
) -> pd.DataFrame:
    max_lag = max(lag_orders)
    ordered = future_df.sort_values(["Date", "series_id"]).copy()
    ordered["row_id"] = np.arange(ordered.shape[0], dtype=int)

    preds: list[tuple[int, float, float]] = []
    for _, block in ordered.groupby("Date", sort=True):
        rows = []
        row_ids = block["row_id"].to_numpy(dtype=int)
        series_ids = block["series_id"].tolist()
        clim_sales = block["clim_sales"].to_numpy(dtype=float)

        for row in block.itertuples(index=False):
            hist = history.get(row.series_id)
            if hist is None:
                hist = deque(maxlen=max_lag)
                history[row.series_id] = hist
            lag_features = {f"lag{lag}": _get_lag_value(hist, lag) for lag in lag_orders}
            rows.append(
                {
                    **lag_features,
                    "week_of_year": float(row.week_of_year),
                    "month": float(row.month),
                    "is_holiday_int": float(row.is_holiday_int),
                    "temp_anom": float(row.temp_anom),
                    "fuel_anom": float(row.fuel_anom),
                    "MarkDown1": float(row.MarkDown1),
                    "MarkDown2": float(row.MarkDown2),
                    "MarkDown3": float(row.MarkDown3),
                    "MarkDown4": float(row.MarkDown4),
                    "MarkDown5": float(row.MarkDown5),
                    "CPI": float(row.CPI),
                    "Unemployment": float(row.Unemployment),
                }
            )

        X_step = pd.DataFrame(rows, columns=feature_cols)
        pred_anom_step = np.asarray(model.predict(X_step), dtype=float)
        pred_anom_step = np.clip(pred_anom_step, lower, upper)
        pred_sales_step = pred_anom_step + clim_sales

        for row_id, series_id, pred_sales, pred_anom in zip(
            row_ids, series_ids, pred_sales_step, pred_anom_step
        ):
            preds.append((int(row_id), float(pred_sales), float(pred_anom)))
            history[series_id].append(float(pred_anom))

    pred_df = pd.DataFrame(preds, columns=["row_id", "pred_sales", "pred_anom"])
    return ordered.merge(pred_df, on="row_id", how="left")

def run_cv_tabular_baseline_anomaly(
    train_df: pd.DataFrame,
    model_builder: Callable[[], object],
    n_folds: int,
    val_weeks: int,
    lag_orders: tuple[int, ...] = LAG_ORDERS,
    feature_cols: list[str] = ANOM_FEATURE_COLS,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    folds = build_forward_folds(train_df, n_folds=n_folds, val_weeks=val_weeks)
    metrics_rows: list[dict[str, float | int | str]] = []
    oof_parts: list[pd.DataFrame] = []
    max_lag = max(lag_orders)

    for fold_id, train_dates, val_dates in folds:
        print(f"\nFold {fold_id}: train={len(train_dates)} weeks, val={len(val_dates)} weeks")
        t0 = time.perf_counter()

        fold_train_raw = train_df[train_df["Date"].isin(train_dates)].copy()
        fold_val_raw = train_df[train_df["Date"].isin(val_dates)].copy()

        climatology = compute_climatology(fold_train_raw)
        fold_train = attach_climatology(fold_train_raw, climatology).sort_values(["series_id", "Date"]).copy()
        fold_val = attach_climatology(fold_val_raw, climatology).sort_values(["Date", "series_id"]).copy()

        X_train, y_train = make_supervised_lag_frame(fold_train, lag_orders=lag_orders, feature_cols=feature_cols)
        model = model_builder()
        model.fit(X_train, y_train)

        history = initialize_target_history(fold_train, max_lag=max_lag)
        lower, upper = _robust_bounds(fold_train[TARGET_COL].to_numpy(dtype=float), clip_scale=2.0)
        fold_val_pred = recursive_predict_tabular_anomaly(
            fold_val,
            model,
            history,
            lag_orders=lag_orders,
            feature_cols=feature_cols,
            lower=lower,
            upper=upper,
        )

        fold_metrics = summarize_metrics(fold_val_pred)
        fold_metrics["runtime_sec"] = float(time.perf_counter() - t0)
        fold_metrics["fold"] = fold_id
        fold_metrics["train_start"] = str(pd.Timestamp(train_dates[0]).date())
        fold_metrics["train_end"] = str(pd.Timestamp(train_dates[-1]).date())
        fold_metrics["val_start"] = str(pd.Timestamp(val_dates[0]).date())
        fold_metrics["val_end"] = str(pd.Timestamp(val_dates[-1]).date())
        metrics_rows.append(fold_metrics)

        print(
            f"  wmae={fold_metrics['wmae']:.4f} mae={fold_metrics['mae']:.4f} "
            f"rmse={fold_metrics['rmse']:.4f} runtime={fold_metrics['runtime_sec']:.1f}s"
        )

        keep_cols = [
            "Store",
            "Dept",
            "Date",
            "IsHoliday",
            "Weekly_Sales",
            "sales_anom",
            "pred_anom",
            "pred_sales",
            "series_id",
        ]
        part = fold_val_pred[keep_cols].copy()
        part["fold"] = fold_id
        oof_parts.append(part)

    metrics_df = pd.DataFrame(metrics_rows)
    if not metrics_df.empty:
        mean_row = {"fold": "mean"}
        for col in ["wmae", "mae", "rmse", "anom_mae", "n_rows", "runtime_sec"]:
            if col in metrics_df.columns:
                mean_row[col] = float(metrics_df[col].mean())
        metrics_df = pd.concat([metrics_df, pd.DataFrame([mean_row])], ignore_index=True)

    oof_df = pd.concat(oof_parts, ignore_index=True) if oof_parts else pd.DataFrame()
    return metrics_df, oof_df

def fit_full_train_and_predict_test_tabular_anomaly(
    train_df: pd.DataFrame,
    test_df: pd.DataFrame,
    model_builder: Callable[[], object],
    lag_orders: tuple[int, ...] = LAG_ORDERS,
    feature_cols: list[str] = ANOM_FEATURE_COLS,
) -> pd.DataFrame:
    climatology = compute_climatology(train_df)
    full_train = attach_climatology(train_df.copy(), climatology).sort_values(["series_id", "Date"]).copy()
    full_test = attach_climatology(test_df.copy(), climatology).sort_values(["Date", "series_id"]).copy()

    X_train, y_train = make_supervised_lag_frame(full_train, lag_orders=lag_orders, feature_cols=feature_cols)
    model = model_builder()
    model.fit(X_train, y_train)

    history = initialize_target_history(full_train, max_lag=max(lag_orders))
    lower, upper = _robust_bounds(full_train[TARGET_COL].to_numpy(dtype=float), clip_scale=2.0)
    test_pred = recursive_predict_tabular_anomaly(
        full_test,
        model,
        history,
        lag_orders=lag_orders,
        feature_cols=feature_cols,
        lower=lower,
        upper=upper,
    )
    return test_pred[
        ["Store", "Dept", "Date", "IsHoliday", "pred_sales", "pred_anom", "series_id"]
    ].copy()

2.4 Phase 2 Core Definitions

import json
import time
from collections import deque
from dataclasses import dataclass
from pathlib import Path

import numpy as np
import pandas as pd
import pymc as pm

TARGET_COL = "sales_anom"
PERIOD_WEEKS = 52.0
LAG_ORDERS = (1, 2)
LAG_COLS = [f"sales_anom_lag{lag}" for lag in LAG_ORDERS]
EXOG_COLS = [
    "temp_anom",
    "fuel_anom",
    "MarkDown1",
    "MarkDown2",
    "MarkDown3",
    "MarkDown4",
    "MarkDown5",
    "CPI",
    "Unemployment",
]

@dataclass(slots=True)
class ScaleStats:
    mean: float
    std: float

@dataclass(slots=True)
class StructuralARConfig:
    n_folds: int
    val_weeks: int
    max_eval: int
    random_seed: int
    pred_draws: int
    fourier_order: int
    sigma_clusters: int

def ensure_dirs(base: Path) -> dict[str, Path]:
    paths = {"root": base, "metrics": base / "metrics", "preds": base / "preds"}
    for p in paths.values():
        p.mkdir(parents=True, exist_ok=True)
    return paths

def fit_standardizer(series: pd.Series) -> ScaleStats:
    mean = float(series.mean())
    std = float(series.std())
    if not np.isfinite(std) or std < 1e-8:
        std = 1.0
    return ScaleStats(mean=mean, std=std)

def zscore(values: np.ndarray, stats: ScaleStats) -> np.ndarray:
    return (values - stats.mean) / stats.std

def add_structural_features(
    train_df: pd.DataFrame,
    val_df: pd.DataFrame,
    fourier_order: int,
) -> tuple[pd.DataFrame, pd.DataFrame, list[str]]:
    min_date = train_df["Date"].min()
    for frame in (train_df, val_df):
        frame["t_week"] = ((frame["Date"] - min_date).dt.days // 7).astype(float)
        frame["trend"] = frame["t_week"]

    fourier_cols: list[str] = []
    for k in range(1, fourier_order + 1):
        sin_col = f"ff_sin_{k}"
        cos_col = f"ff_cos_{k}"
        angle_train = 2.0 * np.pi * k * train_df["t_week"] / PERIOD_WEEKS
        angle_val = 2.0 * np.pi * k * val_df["t_week"] / PERIOD_WEEKS
        train_df[sin_col] = np.sin(angle_train)
        train_df[cos_col] = np.cos(angle_train)
        val_df[sin_col] = np.sin(angle_val)
        val_df[cos_col] = np.cos(angle_val)
        fourier_cols.extend([sin_col, cos_col])

    return train_df, val_df, fourier_cols

def prepare_fold_data(
    train_df: pd.DataFrame,
    train_dates: np.ndarray,
    val_dates: np.ndarray,
    fourier_order: int,
) -> tuple[pd.DataFrame, pd.DataFrame, list[str]]:
    fold_train_raw = train_df[train_df["Date"].isin(train_dates)].copy()
    fold_val_raw = train_df[train_df["Date"].isin(val_dates)].copy()
    fold_train_raw["split"] = "train"
    fold_val_raw["split"] = "val"

    climatology = compute_climatology(fold_train_raw)
    combined = pd.concat([fold_train_raw, fold_val_raw], ignore_index=True)
    combined = attach_climatology(combined, climatology)
    combined = add_anomaly_lags(combined)
    combined = combined.sort_values(["Date", "Store", "Dept"]).reset_index(drop=True)

    fold_train = combined[combined["split"] == "train"].copy()
    fold_val = combined[combined["split"] == "val"].copy()

    # Keep early rows by treating unavailable long-lag anomalies as neutral (0).
    for c in LAG_COLS:
        fold_train[c] = fold_train[c].fillna(0.0)
        fold_val[c] = fold_val[c].fillna(0.0)

    for c in EXOG_COLS:
        median = float(fold_train[c].median())
        if not np.isfinite(median):
            median = 0.0
        fold_train[c] = fold_train[c].fillna(median)
        fold_val[c] = fold_val[c].fillna(median)

    fold_train = fold_train.dropna(subset=[TARGET_COL]).reset_index(drop=True)
    fold_val = fold_val.dropna(subset=[TARGET_COL]).reset_index(drop=True)

    fold_train, fold_val, fourier_cols = add_structural_features(
        fold_train, fold_val, fourier_order=fourier_order
    )
    return fold_train, fold_val, fourier_cols

def _lag_from_history(hist: deque[float], lag: int) -> float:
    if len(hist) >= lag:
        return float(hist[-lag])
    return 0.0

def initialize_history(train_part: pd.DataFrame) -> dict[str, deque[float]]:
    history: dict[str, deque[float]] = {}
    max_lag = max(LAG_ORDERS)
    ordered = train_part.sort_values(["series_id", "Date"])
    for series_id, grp in ordered.groupby("series_id", sort=False):
        vals = grp[TARGET_COL].dropna().to_numpy(dtype=float)
        history[series_id] = deque(vals.tolist(), maxlen=max_lag)
    return history

def build_sigma_clusters(train_part: pd.DataFrame, n_clusters: int) -> dict[str, int]:
    if n_clusters <= 1:
        return {}
    series_stats = (
        train_part.groupby("series_id", as_index=False)[TARGET_COL]
        .std()
        .rename(columns={TARGET_COL: "anom_std"})
    )
    series_stats["anom_std"] = series_stats["anom_std"].fillna(0.0)
    k = int(max(1, min(n_clusters, series_stats.shape[0])))
    if k == 1:
        return {sid: 0 for sid in series_stats["series_id"].tolist()}
    # Quantile bins by anomaly volatility to capture heteroskedastic groups.
    labels = pd.qcut(series_stats["anom_std"], q=k, labels=False, duplicates="drop")
    labels = labels.fillna(0).astype(int)
    return dict(zip(series_stats["series_id"], labels))

def fit_structural_model(
    train_part: pd.DataFrame,
    fourier_cols: list[str],
    cfg: StructuralARConfig,
) -> dict:
    y_stats = fit_standardizer(train_part[TARGET_COL])
    lag_stats = [fit_standardizer(train_part[c]) for c in LAG_COLS]
    exog_stats = [fit_standardizer(train_part[c]) for c in EXOG_COLS]
    trend_stats = fit_standardizer(train_part["trend"])

    y_train = zscore(train_part[TARGET_COL].to_numpy(dtype=float), y_stats)
    lag_z_matrix = np.column_stack(
        [
            zscore(train_part[col].to_numpy(dtype=float), stats)
            for col, stats in zip(LAG_COLS, lag_stats)
        ]
    )
    exog_z_matrix = np.column_stack(
        [
            zscore(train_part[col].to_numpy(dtype=float), stats)
            for col, stats in zip(EXOG_COLS, exog_stats)
        ]
    )
    trend_z = zscore(train_part["trend"].to_numpy(dtype=float), trend_stats)
    holiday = train_part["is_holiday_int"].to_numpy(dtype=float)
    ff_matrix = train_part[fourier_cols].to_numpy(dtype=float)

    series_ids = train_part["series_id"].drop_duplicates().sort_values().tolist()
    series_to_idx = {sid: i for i, sid in enumerate(series_ids)}
    train_series_idx = train_part["series_id"].map(series_to_idx).to_numpy(dtype=int)
    n_series = len(series_ids)
    series_to_cluster: dict[str, int] = {}
    cluster_idx: np.ndarray | None = None
    n_cluster = 0
    if cfg.sigma_clusters > 1:
        series_to_cluster = build_sigma_clusters(train_part, cfg.sigma_clusters)
        if series_to_cluster:
            cluster_idx = train_part["series_id"].map(series_to_cluster).fillna(0).to_numpy(dtype=int)
            n_cluster = int(np.max(cluster_idx)) + 1

    with pm.Model() as model:
        alpha_mu = pm.Normal("alpha_mu", mu=0.0, sigma=1.0)
        alpha_sigma = pm.HalfNormal("alpha_sigma", sigma=1.0)
        z_alpha = pm.Normal("z_alpha", mu=0.0, sigma=1.0, shape=n_series)
        alpha_series = pm.Deterministic("alpha_series", alpha_mu + alpha_sigma * z_alpha)

        beta_lags = pm.Normal("beta_lags", mu=0.0, sigma=0.6, shape=len(LAG_COLS))
        beta_exog = pm.Normal("beta_exog", mu=0.0, sigma=0.5, shape=len(EXOG_COLS))
        beta_holiday = pm.Normal("beta_holiday", mu=0.0, sigma=0.8)
        beta_trend = pm.Normal("beta_trend", mu=0.0, sigma=0.5)
        beta_fourier = pm.Normal("beta_fourier", mu=0.0, sigma=0.5, shape=ff_matrix.shape[1])
        if cluster_idx is not None and n_cluster > 1:
            sigma_cluster = pm.HalfNormal("sigma_cluster", sigma=0.8, shape=n_cluster)
            sigma_obs = sigma_cluster[cluster_idx]
        else:
            sigma = pm.HalfNormal("sigma", sigma=0.8)
            sigma_obs = sigma

        mu = (
            alpha_series[train_series_idx]
            + pm.math.dot(lag_z_matrix, beta_lags)
            + pm.math.dot(exog_z_matrix, beta_exog)
            + beta_holiday * holiday
            + beta_trend * trend_z
            + pm.math.dot(ff_matrix, beta_fourier)
        )
        pm.Normal("y_obs", mu=mu, sigma=sigma_obs, observed=y_train)

        map_point = pm.find_MAP(
            progressbar=False,
            maxeval=cfg.max_eval,
            return_raw=False,
            seed=cfg.random_seed,
        )

    train_vals = train_part[TARGET_COL].to_numpy(dtype=float)
    q01, q99 = np.quantile(train_vals, [0.01, 0.99])
    span = max(float(q99 - q01), float(np.std(train_vals)), 1.0)
    lower = float(q01 - 2.0 * span)
    upper = float(q99 + 2.0 * span)

    fitted = {
        "map_point": map_point,
        "series_to_idx": series_to_idx,
        "y_stats": y_stats,
        "lag_stats": lag_stats,
        "exog_stats": exog_stats,
        "exog_cols": EXOG_COLS,
        "trend_stats": trend_stats,
        "fourier_cols": fourier_cols,
        "history_base": initialize_history(train_part),
        "clip_bounds": (lower, upper),
        "series_to_cluster": series_to_cluster,
    }
    return fitted

def recursive_predict_with_uncertainty(
    val_part: pd.DataFrame,
    fitted: dict,
    cfg: StructuralARConfig,
) -> tuple[np.ndarray, np.ndarray]:
    ordered = (
        val_part.reset_index()
        .rename(columns={"index": "orig_idx"})
        .sort_values(["Date", "series_id"])
        .reset_index(drop=True)
    )

    mp = fitted["map_point"]
    series_to_idx = fitted["series_to_idx"]
    alpha_series = np.asarray(mp["alpha_series"], dtype=float)
    alpha_mu = float(mp["alpha_mu"])

    beta_lags = np.asarray(mp["beta_lags"], dtype=float)
    beta_exog = np.asarray(mp["beta_exog"], dtype=float)
    beta_holiday = float(mp["beta_holiday"])
    beta_trend = float(mp["beta_trend"])
    beta_fourier = np.asarray(mp["beta_fourier"], dtype=float)
    sigma_global = float(mp["sigma"]) if "sigma" in mp else None
    sigma_cluster = np.asarray(mp["sigma_cluster"], dtype=float) if "sigma_cluster" in mp else None

    y_stats = fitted["y_stats"]
    lag_stats = fitted["lag_stats"]
    exog_stats = fitted["exog_stats"]
    exog_cols = fitted["exog_cols"]
    trend_stats = fitted["trend_stats"]
    fourier_cols = fitted["fourier_cols"]
    series_to_cluster = fitted.get("series_to_cluster", {})
    lower, upper = fitted["clip_bounds"]

    sid_arr = ordered["series_id"].to_numpy()
    exog_arr = ordered[exog_cols].to_numpy(dtype=float)
    holiday_arr = ordered["is_holiday_int"].to_numpy(dtype=float)
    trend_arr = ordered["trend"].to_numpy(dtype=float)
    ff_arr = ordered[fourier_cols].to_numpy(dtype=float)
    orig_idx = ordered["orig_idx"].to_numpy(dtype=int)

    n = len(ordered)
    draws = max(1, cfg.pred_draws)
    draw_preds = np.empty((draws, n), dtype=float)
    rng = np.random.default_rng(cfg.random_seed)

    for d in range(draws):
        history = {sid: deque(vals, maxlen=max(LAG_ORDERS)) for sid, vals in fitted["history_base"].items()}
        for i in range(n):
            sid = sid_arr[i]
            hist = history.get(sid)
            if hist is None:
                hist = deque(maxlen=max(LAG_ORDERS))
                history[sid] = hist
            sid_idx = series_to_idx.get(sid)
            alpha = alpha_series[sid_idx] if sid_idx is not None else alpha_mu

            lag_z_vals = [
                (_lag_from_history(hist, lag) - stats.mean) / stats.std
                for lag, stats in zip(LAG_ORDERS, lag_stats)
            ]
            exog_z_vals = [
                (exog_arr[i, j] - stats.mean) / stats.std
                for j, stats in enumerate(exog_stats)
            ]
            trend_z = (trend_arr[i] - trend_stats.mean) / trend_stats.std

            pred_z = (
                alpha
                + float(np.dot(lag_z_vals, beta_lags))
                + float(np.dot(exog_z_vals, beta_exog))
                + beta_holiday * holiday_arr[i]
                + beta_trend * trend_z
                + float(np.dot(ff_arr[i], beta_fourier))
            )
            if sigma_cluster is not None and series_to_cluster:
                c_idx = int(series_to_cluster.get(sid, 0))
                sigma_draw = float(sigma_cluster[min(c_idx, len(sigma_cluster) - 1)])
            else:
                sigma_draw = float(sigma_global if sigma_global is not None else 1.0)
            pred_z += float(rng.normal(loc=0.0, scale=sigma_draw))
            pred_raw = float(pred_z * y_stats.std + y_stats.mean)
            pred_raw = float(np.clip(pred_raw, lower, upper))

            draw_preds[d, i] = pred_raw
            hist.append(pred_raw)

    mean_ordered = draw_preds.mean(axis=0)
    sd_ordered = draw_preds.std(axis=0, ddof=0)
    mean_series = pd.Series(mean_ordered, index=orig_idx)
    sd_series = pd.Series(sd_ordered, index=orig_idx)
    pred_mean = mean_series.reindex(val_part.index).to_numpy(dtype=float)
    pred_sd = sd_series.reindex(val_part.index).to_numpy(dtype=float)
    return pred_mean, pred_sd

def evaluate_fold(val_df: pd.DataFrame) -> dict[str, float]:
    y_true = val_df["Weekly_Sales"].to_numpy(dtype=float)
    y_pred = val_df["pred_sales"].to_numpy(dtype=float)
    holiday = val_df["IsHoliday"].to_numpy(dtype=bool)
    return {
        "wmae": wmae(y_true, y_pred, holiday),
        "mae": float(np.mean(np.abs(y_true - y_pred))),
        "rmse": float(np.sqrt(np.mean((y_true - y_pred) ** 2))),
        "anom_mae": float(np.mean(np.abs(val_df["sales_anom"] - val_df["pred_anom"]))),
        "n_rows": float(len(val_df)),
    }

def run_cv(train_df: pd.DataFrame, cfg: StructuralARConfig) -> tuple[pd.DataFrame, pd.DataFrame]:
    folds = build_forward_folds(train_df, n_folds=cfg.n_folds, val_weeks=cfg.val_weeks)
    metrics_rows: list[dict] = []
    oof_parts: list[pd.DataFrame] = []

    for fold_id, train_dates, val_dates in folds:
        print(f"\nFold {fold_id}: train={len(train_dates)} weeks, val={len(val_dates)} weeks")
        fold_train, fold_val, fourier_cols = prepare_fold_data(
            train_df, train_dates, val_dates, fourier_order=cfg.fourier_order
        )
        print(f"  rows used -> train={len(fold_train)}, val={len(fold_val)}")

        t0 = time.perf_counter()
        fitted = fit_structural_model(fold_train, fourier_cols=fourier_cols, cfg=cfg)
        pred_anom, pred_anom_sd = recursive_predict_with_uncertainty(fold_val, fitted=fitted, cfg=cfg)
        elapsed = time.perf_counter() - t0

        fold_val = fold_val.copy()
        fold_val["pred_anom"] = pred_anom
        fold_val["pred_anom_sd"] = pred_anom_sd
        fold_val["pred_sales"] = fold_val["pred_anom"] + fold_val["clim_sales"]
        fold_val["pred_sales_sd"] = fold_val["pred_anom_sd"]
        fold_val["pred_sales_lower"] = fold_val["pred_sales"] - 1.96 * fold_val["pred_sales_sd"]
        fold_val["pred_sales_upper"] = fold_val["pred_sales"] + 1.96 * fold_val["pred_sales_sd"]
        fold_val["fold"] = fold_id

        fold_metrics = evaluate_fold(fold_val)
        fold_metrics["fold"] = fold_id
        fold_metrics["train_start"] = str(pd.Timestamp(train_dates[0]).date())
        fold_metrics["train_end"] = str(pd.Timestamp(train_dates[-1]).date())
        fold_metrics["val_start"] = str(pd.Timestamp(val_dates[0]).date())
        fold_metrics["val_end"] = str(pd.Timestamp(val_dates[-1]).date())
        fold_metrics["runtime_sec"] = float(elapsed)
        metrics_rows.append(fold_metrics)
        print(
            f"  wmae={fold_metrics['wmae']:.4f} mae={fold_metrics['mae']:.4f} "
            f"rmse={fold_metrics['rmse']:.4f} runtime={elapsed:.1f}s"
        )

        keep_cols = [
            "Store",
            "Dept",
            "Date",
            "IsHoliday",
            "Weekly_Sales",
            "sales_anom",
            "pred_anom",
            "pred_anom_sd",
            "pred_sales",
            "pred_sales_sd",
            "pred_sales_lower",
            "pred_sales_upper",
            "fold",
        ]
        oof_parts.append(fold_val[keep_cols].copy())

    metrics_df = pd.DataFrame(metrics_rows)
    if not metrics_df.empty:
        mean_row = {"fold": "mean"}
        for col in ["wmae", "mae", "rmse", "anom_mae", "n_rows", "runtime_sec"]:
            mean_row[col] = float(metrics_df[col].mean())
        metrics_df = pd.concat([metrics_df, pd.DataFrame([mean_row])], ignore_index=True)

    oof_df = pd.concat(oof_parts, ignore_index=True) if oof_parts else pd.DataFrame()
    return metrics_df, oof_df


def phase2_default_config() -> dict:
    return {
        "train_path": "train_feat.parquet",
        "test_path": "test_feat.parquet",
        "output_dir": "outputs/phase2_structural_ar_full",
        "n_folds": 4,
        "val_weeks": 13,
        "max_eval": 1800,
        "random_seed": 8927,
        "pred_draws": 40,
        "fourier_order": 3,
        "sigma_clusters": 0,
        "max_series": 3331,
    }


def run_phase2(
    config: dict | None = None,
    _ensure_dirs=ensure_dirs,
    _StructuralARConfig=StructuralARConfig,
    _load_data=load_data,
    _run_cv=run_cv,
) -> dict[str, Path]:
    cfg = phase2_default_config()
    if config:
        cfg.update(config)

    out = _ensure_dirs(Path(cfg["output_dir"]))
    model_cfg = _StructuralARConfig(
        n_folds=int(cfg["n_folds"]),
        val_weeks=int(cfg["val_weeks"]),
        max_eval=int(cfg["max_eval"]),
        random_seed=int(cfg["random_seed"]),
        pred_draws=int(cfg["pred_draws"]),
        fourier_order=int(cfg["fourier_order"]),
        sigma_clusters=int(cfg["sigma_clusters"]),
    )

    max_series = cfg.get("max_series")
    if max_series is not None:
        max_series = int(max_series)

    train_df, _ = _load_data(
        Path(cfg["train_path"]),
        Path(cfg["test_path"]),
        max_series=max_series,
    )

    metrics_df, oof_df = _run_cv(train_df, cfg=model_cfg)
    metrics_path = out["metrics"] / "structural_ar_cv.csv"
    oof_path = out["preds"] / "structural_ar_oof.parquet"
    metrics_df.to_csv(metrics_path, index=False)
    oof_df.to_parquet(oof_path, index=False)

    config_payload = {
        **cfg,
        "lag_orders": list(LAG_ORDERS),
        "exogenous_features": EXOG_COLS,
        "note": "Structural AR model with hierarchical series intercept, trend, Fourier seasonality, and recursive validation.",
    }
    config_path = out["metrics"] / "structural_ar_config.json"
    config_path.write_text(json.dumps(config_payload, indent=2), encoding="utf-8")

    return {
        "metrics": metrics_path,
        "oof": oof_path,
        "config": config_path,
    }


def tune_phase2(
    base_config: dict,
    search_space: list[dict],
    tune_overrides: dict | None = None,
    search_output_dir: str = "outputs/phase2_structural_ar_full/search",
    _StructuralARConfig=StructuralARConfig,
    _load_data=load_data,
    _run_cv=run_cv,
) -> tuple[pd.DataFrame, dict]:
    cfg0 = phase2_default_config()
    cfg0.update(base_config)
    if tune_overrides:
        cfg0.update(tune_overrides)

    rows: list[dict] = []
    best_cfg: dict | None = None
    best_wmae = float("inf")

    for idx, cand in enumerate(search_space, start=1):
        cfg = {**cfg0, **cand}
        model_cfg = _StructuralARConfig(
            n_folds=int(cfg["n_folds"]),
            val_weeks=int(cfg["val_weeks"]),
            max_eval=int(cfg["max_eval"]),
            random_seed=int(cfg["random_seed"]),
            pred_draws=int(cfg["pred_draws"]),
            fourier_order=int(cfg["fourier_order"]),
            sigma_clusters=int(cfg["sigma_clusters"]),
        )

        max_series = cfg.get("max_series")
        if max_series is not None:
            max_series = int(max_series)

        train_df, _ = _load_data(
            Path(cfg["train_path"]),
            Path(cfg["test_path"]),
            max_series=max_series,
        )

        metrics_df, _ = _run_cv(train_df, cfg=model_cfg)
        mean_row = mean_row_from_metrics(metrics_df)
        wmae_val = float(mean_row["wmae"])
        rec = {
            "candidate": idx,
            **cand,
            "wmae": wmae_val,
            "mae": float(mean_row.get("mae", np.nan)),
            "rmse": float(mean_row.get("rmse", np.nan)),
        }
        rows.append(rec)
        if wmae_val < best_wmae:
            best_wmae = wmae_val
            best_cfg = cfg

    results_df = pd.DataFrame(rows).sort_values("wmae").reset_index(drop=True)
    out_dir = Path(search_output_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    results_df.to_csv(out_dir / "search_results.csv", index=False)
    if best_cfg is None:
        raise RuntimeError("Phase 2 tuning did not evaluate any candidates")
    (out_dir / "best_config.json").write_text(json.dumps(best_cfg, indent=2), encoding="utf-8")
    return results_df, best_cfg

2.5 Phase 3 Core Definitions

import json
from pathlib import Path

from sklearn.impute import SimpleImputer
from sklearn.linear_model import ElasticNet
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

def make_model(alpha: float, l1_ratio: float, max_iter: int):
    return Pipeline(
        steps=[
            ("imputer", SimpleImputer(strategy="median")),
            ("scaler", StandardScaler()),
            (
                "model",
                ElasticNet(
                    alpha=alpha,
                    l1_ratio=l1_ratio,
                    max_iter=max_iter,
                    random_state=42,
                ),
            ),
        ]
    )


def phase3_default_config() -> dict:
    return {
        "train_path": "train_feat.parquet",
        "test_path": "test_feat.parquet",
        "output_dir": "outputs/baselines/elastic_net",
        "n_folds": 4,
        "val_weeks": 13,
        "max_series": None,
        "alpha": 0.02,
        "l1_ratio": 0.2,
        "max_iter": 10000,
    }


def run_phase3(
    config: dict | None = None,
    _ensure_output_dir=ensure_output_dir,
    _load_data=load_data,
    _run_cv=run_cv_tabular_baseline_anomaly,
    _fit_full=fit_full_train_and_predict_test_tabular_anomaly,
    _make_model=make_model,
) -> dict[str, Path]:
    cfg = phase3_default_config()
    if config:
        cfg.update(config)

    out_dir = _ensure_output_dir(Path(cfg["output_dir"]))
    train_df, test_df = _load_data(
        Path(cfg["train_path"]),
        Path(cfg["test_path"]),
        max_series=cfg.get("max_series"),
    )

    metrics_df, oof_df = _run_cv(
        train_df,
        model_builder=lambda: _make_model(float(cfg["alpha"]), float(cfg["l1_ratio"]), int(cfg["max_iter"])),
        n_folds=int(cfg["n_folds"]),
        val_weeks=int(cfg["val_weeks"]),
        lag_orders=LAG_ORDERS,
        feature_cols=ANOM_FEATURE_COLS,
    )

    metrics_path = out_dir / "metrics.csv"
    oof_path = out_dir / "oof.parquet"
    metrics_df.to_csv(metrics_path, index=False)
    oof_df.to_parquet(oof_path, index=False)

    test_pred = _fit_full(
        train_df,
        test_df,
        model_builder=lambda: _make_model(float(cfg["alpha"]), float(cfg["l1_ratio"]), int(cfg["max_iter"])),
        lag_orders=LAG_ORDERS,
        feature_cols=ANOM_FEATURE_COLS,
    )
    test_path = out_dir / "test.parquet"
    test_pred.to_parquet(test_path, index=False)

    config_payload = {
        **cfg,
        "model": "elastic_net",
        "target": "sales_anom",
        "features": ANOM_FEATURE_COLS,
        "lags": list(LAG_ORDERS),
    }
    config_path = out_dir / "config.json"
    config_path.write_text(json.dumps(config_payload, indent=2), encoding="utf-8")

    return {"metrics": metrics_path, "oof": oof_path, "test": test_path, "config": config_path}


def tune_phase3(
    base_config: dict,
    search_space: list[dict],
    search_output_dir: str = "outputs/baselines/elastic_net/search",
    _load_data=load_data,
    _run_cv=run_cv_tabular_baseline_anomaly,
    _make_model=make_model,
) -> tuple[pd.DataFrame, dict]:
    cfg0 = phase3_default_config()
    cfg0.update(base_config)
    train_df, _ = _load_data(Path(cfg0["train_path"]), Path(cfg0["test_path"]), max_series=cfg0.get("max_series"))

    rows: list[dict] = []
    best_cfg: dict | None = None
    best_wmae = float("inf")

    for idx, cand in enumerate(search_space, start=1):
        cfg = {**cfg0, **cand}
        metrics_df, _ = _run_cv(
            train_df,
            model_builder=lambda: _make_model(float(cfg["alpha"]), float(cfg["l1_ratio"]), int(cfg["max_iter"])),
            n_folds=int(cfg["n_folds"]),
            val_weeks=int(cfg["val_weeks"]),
            lag_orders=LAG_ORDERS,
            feature_cols=ANOM_FEATURE_COLS,
        )
        mean_row = mean_row_from_metrics(metrics_df)
        wmae_val = float(mean_row["wmae"])
        rows.append(
            {
                "candidate": idx,
                **cand,
                "wmae": wmae_val,
                "mae": float(mean_row.get("mae", np.nan)),
                "rmse": float(mean_row.get("rmse", np.nan)),
            }
        )
        if wmae_val < best_wmae:
            best_wmae = wmae_val
            best_cfg = cfg

    results_df = pd.DataFrame(rows).sort_values("wmae").reset_index(drop=True)
    out_dir = Path(search_output_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    results_df.to_csv(out_dir / "search_results.csv", index=False)
    if best_cfg is None:
        raise RuntimeError("Phase 3 tuning did not evaluate any candidates")
    (out_dir / "best_config.json").write_text(json.dumps(best_cfg, indent=2), encoding="utf-8")
    return results_df, best_cfg

2.6 Phase 4 Core Definitions

import json
from pathlib import Path

from sklearn.ensemble import RandomForestRegressor
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline

def make_model(n_estimators: int, max_depth: int, min_samples_leaf: int, max_features: str):
    return Pipeline(
        steps=[
            ("imputer", SimpleImputer(strategy="median")),
            (
                "model",
                RandomForestRegressor(
                    n_estimators=n_estimators,
                    max_depth=max_depth,
                    min_samples_leaf=min_samples_leaf,
                    max_features=max_features,
                    random_state=42,
                    n_jobs=-1,
                ),
            ),
        ]
    )


def phase4_default_config() -> dict:
    return {
        "train_path": "train_feat.parquet",
        "test_path": "test_feat.parquet",
        "output_dir": "outputs/baselines/random_forest",
        "n_folds": 4,
        "val_weeks": 13,
        "max_series": None,
        "n_estimators": 120,
        "max_depth": 18,
        "min_samples_leaf": 2,
        "max_features": "sqrt",
    }


def run_phase4(
    config: dict | None = None,
    _ensure_output_dir=ensure_output_dir,
    _load_data=load_data,
    _run_cv=run_cv_tabular_baseline_anomaly,
    _fit_full=fit_full_train_and_predict_test_tabular_anomaly,
    _make_model=make_model,
) -> dict[str, Path]:
    cfg = phase4_default_config()
    if config:
        cfg.update(config)

    out_dir = _ensure_output_dir(Path(cfg["output_dir"]))
    train_df, test_df = _load_data(Path(cfg["train_path"]), Path(cfg["test_path"]), max_series=cfg.get("max_series"))

    metrics_df, oof_df = _run_cv(
        train_df,
        model_builder=lambda: _make_model(
            int(cfg["n_estimators"]),
            int(cfg["max_depth"]),
            int(cfg["min_samples_leaf"]),
            str(cfg["max_features"]),
        ),
        n_folds=int(cfg["n_folds"]),
        val_weeks=int(cfg["val_weeks"]),
        lag_orders=LAG_ORDERS,
        feature_cols=ANOM_FEATURE_COLS,
    )

    metrics_path = out_dir / "metrics.csv"
    oof_path = out_dir / "oof.parquet"
    metrics_df.to_csv(metrics_path, index=False)
    oof_df.to_parquet(oof_path, index=False)

    test_pred = _fit_full(
        train_df,
        test_df,
        model_builder=lambda: _make_model(
            int(cfg["n_estimators"]),
            int(cfg["max_depth"]),
            int(cfg["min_samples_leaf"]),
            str(cfg["max_features"]),
        ),
        lag_orders=LAG_ORDERS,
        feature_cols=ANOM_FEATURE_COLS,
    )
    test_path = out_dir / "test.parquet"
    test_pred.to_parquet(test_path, index=False)

    config_payload = {
        **cfg,
        "model": "random_forest",
        "target": "sales_anom",
        "features": ANOM_FEATURE_COLS,
        "lags": list(LAG_ORDERS),
    }
    config_path = out_dir / "config.json"
    config_path.write_text(json.dumps(config_payload, indent=2), encoding="utf-8")

    return {"metrics": metrics_path, "oof": oof_path, "test": test_path, "config": config_path}


def tune_phase4(
    base_config: dict,
    search_space: list[dict],
    search_output_dir: str = "outputs/baselines/random_forest/search",
    _load_data=load_data,
    _run_cv=run_cv_tabular_baseline_anomaly,
    _make_model=make_model,
) -> tuple[pd.DataFrame, dict]:
    cfg0 = phase4_default_config()
    cfg0.update(base_config)
    train_df, _ = _load_data(Path(cfg0["train_path"]), Path(cfg0["test_path"]), max_series=cfg0.get("max_series"))

    rows: list[dict] = []
    best_cfg: dict | None = None
    best_wmae = float("inf")

    for idx, cand in enumerate(search_space, start=1):
        cfg = {**cfg0, **cand}
        metrics_df, _ = _run_cv(
            train_df,
            model_builder=lambda: _make_model(
                int(cfg["n_estimators"]),
                int(cfg["max_depth"]),
                int(cfg["min_samples_leaf"]),
                str(cfg["max_features"]),
            ),
            n_folds=int(cfg["n_folds"]),
            val_weeks=int(cfg["val_weeks"]),
            lag_orders=LAG_ORDERS,
            feature_cols=ANOM_FEATURE_COLS,
        )
        mean_row = mean_row_from_metrics(metrics_df)
        wmae_val = float(mean_row["wmae"])
        rows.append(
            {
                "candidate": idx,
                **cand,
                "wmae": wmae_val,
                "mae": float(mean_row.get("mae", np.nan)),
                "rmse": float(mean_row.get("rmse", np.nan)),
            }
        )
        if wmae_val < best_wmae:
            best_wmae = wmae_val
            best_cfg = cfg

    results_df = pd.DataFrame(rows).sort_values("wmae").reset_index(drop=True)
    out_dir = Path(search_output_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    results_df.to_csv(out_dir / "search_results.csv", index=False)
    if best_cfg is None:
        raise RuntimeError("Phase 4 tuning did not evaluate any candidates")
    (out_dir / "best_config.json").write_text(json.dumps(best_cfg, indent=2), encoding="utf-8")
    return results_df, best_cfg

2.7 Phase 5 Core Definitions

import json
import time
import warnings
from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.impute import SimpleImputer
from sklearn.linear_model import Ridge
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from statsmodels.tools.sm_exceptions import ConvergenceWarning
from statsmodels.tsa.holtwinters import ExponentialSmoothing

EXOG_FEATURE_COLS = [
    "temp_anom",
    "fuel_anom",
    "MarkDown1",
    "MarkDown2",
    "MarkDown3",
    "MarkDown4",
    "MarkDown5",
    "CPI",
    "Unemployment",
]

def make_exog_model(alpha: float) -> Pipeline:
    return Pipeline(
        steps=[
            ("imputer", SimpleImputer(strategy="median")),
            ("scaler", StandardScaler()),
            ("model", Ridge(alpha=alpha)),
        ]
    )

def _safe_score(aic: float | None, sse: float | None) -> float:
    if aic is not None and np.isfinite(aic):
        return float(aic)
    if sse is not None and np.isfinite(sse):
        return float(sse)
    return float("inf")

def forecast_ets_with_fallback(
    train_values: np.ndarray,
    horizon: int,
    seasonal_periods: int,
    fallback_value: float,
) -> tuple[np.ndarray, str]:
    y = np.asarray(train_values, dtype=float)
    y = y[np.isfinite(y)]

    if horizon <= 0:
        return np.empty(0, dtype=float), "empty-horizon"

    if y.size == 0:
        return np.full(horizon, fallback_value, dtype=float), "global-fallback"

    if y.size < 3:
        return np.full(horizon, float(y[-1]), dtype=float), "last-value-fallback"

    candidate_configs: list[dict[str, object]] = [
        {"trend": None, "seasonal": None, "damped_trend": False, "tag": "level"},
        {"trend": "add", "seasonal": None, "damped_trend": False, "tag": "level-trend"},
    ]
    if y.size >= seasonal_periods + 4:
        candidate_configs.append(
            {"trend": None, "seasonal": "add", "damped_trend": False, "tag": "level-seasonal"}
        )
        candidate_configs.append(
            {"trend": "add", "seasonal": "add", "damped_trend": True, "tag": "trend-seasonal-damped"}
        )

    best_forecast: np.ndarray | None = None
    best_score = float("inf")
    best_tag = "last-value-fallback"

    for cfg in candidate_configs:
        try:
            model = ExponentialSmoothing(
                y,
                trend=cfg["trend"],
                seasonal=cfg["seasonal"],
                damped_trend=bool(cfg["damped_trend"]),
                seasonal_periods=seasonal_periods if cfg["seasonal"] else None,
                initialization_method="estimated",
            )
            fit = model.fit(optimized=True, use_brute=False)
            pred = np.asarray(fit.forecast(horizon), dtype=float)
            score = _safe_score(getattr(fit, "aic", None), getattr(fit, "sse", None))
            if np.all(np.isfinite(pred)) and score < best_score:
                best_score = score
                best_forecast = pred
                best_tag = str(cfg["tag"])
        except Exception:
            continue

    if best_forecast is None:
        return np.full(horizon, float(y[-1]), dtype=float), "last-value-fallback"
    return best_forecast, best_tag

def run_cv(
    train_df: pd.DataFrame,
    n_folds: int,
    val_weeks: int,
    seasonal_periods: int,
    exog_alpha: float,
) -> tuple[pd.DataFrame, pd.DataFrame, dict[str, int]]:
    folds = build_forward_folds(train_df, n_folds=n_folds, val_weeks=val_weeks)
    metrics_rows: list[dict[str, float | int | str]] = []
    oof_parts: list[pd.DataFrame] = []
    mode_counts: dict[str, int] = {}

    for fold_id, train_dates, val_dates in folds:
        print(f"\nFold {fold_id}: train={len(train_dates)} weeks, val={len(val_dates)} weeks")
        t0 = time.perf_counter()

        fold_train_raw = train_df[train_df["Date"].isin(train_dates)].copy()
        fold_val_raw = train_df[train_df["Date"].isin(val_dates)].copy()
        climatology = compute_climatology(fold_train_raw)
        fold_train = attach_climatology(fold_train_raw, climatology)
        fold_val = attach_climatology(fold_val_raw, climatology)
        exog_model = make_exog_model(exog_alpha)
        exog_model.fit(fold_train[EXOG_FEATURE_COLS], fold_train["sales_anom"].to_numpy(dtype=float))
        fold_train["exog_anom_pred"] = np.asarray(exog_model.predict(fold_train[EXOG_FEATURE_COLS]), dtype=float)
        fold_val["exog_anom_pred"] = np.asarray(exog_model.predict(fold_val[EXOG_FEATURE_COLS]), dtype=float)
        fold_train["resid_anom"] = fold_train["sales_anom"] - fold_train["exog_anom_pred"]
        fallback = float(fold_train["resid_anom"].median())

        fold_preds: list[pd.DataFrame] = []
        val_series_ids = fold_val["series_id"].drop_duplicates().tolist()
        total_series = len(val_series_ids)
        for i, series_id in enumerate(val_series_ids, start=1):
            train_series = fold_train[fold_train["series_id"] == series_id].sort_values("Date")
            val_series = fold_val[fold_val["series_id"] == series_id].sort_values("Date").copy()
            horizon = val_series.shape[0]

            if train_series.empty:
                pred = np.full(horizon, fallback, dtype=float)
                mode = "global-fallback"
            else:
                pred, mode = forecast_ets_with_fallback(
                    train_series["resid_anom"].to_numpy(dtype=float),
                    horizon=horizon,
                    seasonal_periods=seasonal_periods,
                    fallback_value=fallback,
                )
            mode_counts[mode] = mode_counts.get(mode, 0) + 1
            val_series["pred_resid_anom"] = pred
            val_series["pred_anom"] = val_series["exog_anom_pred"].to_numpy(dtype=float) + val_series[
                "pred_resid_anom"
            ].to_numpy(dtype=float)
            val_series["pred_sales"] = val_series["pred_anom"] + val_series["clim_sales"]
            fold_preds.append(val_series)

            if i % 400 == 0:
                print(f"  processed {i}/{total_series} series...")

        fold_val_pred = pd.concat(fold_preds, ignore_index=True) if fold_preds else fold_val.assign(pred_sales=np.nan)
        fold_metrics = summarize_metrics(fold_val_pred)
        fold_metrics["runtime_sec"] = float(time.perf_counter() - t0)
        fold_metrics["fold"] = fold_id
        fold_metrics["train_start"] = str(pd.Timestamp(train_dates[0]).date())
        fold_metrics["train_end"] = str(pd.Timestamp(train_dates[-1]).date())
        fold_metrics["val_start"] = str(pd.Timestamp(val_dates[0]).date())
        fold_metrics["val_end"] = str(pd.Timestamp(val_dates[-1]).date())
        metrics_rows.append(fold_metrics)

        print(
            f"  wmae={fold_metrics['wmae']:.4f} mae={fold_metrics['mae']:.4f} "
            f"rmse={fold_metrics['rmse']:.4f} runtime={fold_metrics['runtime_sec']:.1f}s"
        )

        keep_cols = [
            "Store",
            "Dept",
            "Date",
            "IsHoliday",
            "Weekly_Sales",
            "sales_anom",
            "pred_anom",
            "pred_sales",
            "series_id",
        ]
        part = fold_val_pred[keep_cols].copy()
        part["fold"] = fold_id
        oof_parts.append(part)

    metrics_df = pd.DataFrame(metrics_rows)
    if not metrics_df.empty:
        mean_row = {"fold": "mean"}
        for col in ["wmae", "mae", "rmse", "anom_mae", "n_rows", "runtime_sec"]:
            mean_row[col] = float(metrics_df[col].mean())
        metrics_df = pd.concat([metrics_df, pd.DataFrame([mean_row])], ignore_index=True)

    oof_df = pd.concat(oof_parts, ignore_index=True) if oof_parts else pd.DataFrame()
    return metrics_df, oof_df, mode_counts

def fit_full_and_predict_test(
    train_df: pd.DataFrame,
    test_df: pd.DataFrame,
    seasonal_periods: int,
    exog_alpha: float,
) -> tuple[pd.DataFrame, dict[str, int]]:
    climatology = compute_climatology(train_df)
    full_train = attach_climatology(train_df.copy(), climatology)
    full_test = attach_climatology(test_df.copy(), climatology)
    exog_model = make_exog_model(exog_alpha)
    exog_model.fit(full_train[EXOG_FEATURE_COLS], full_train["sales_anom"].to_numpy(dtype=float))
    full_train["exog_anom_pred"] = np.asarray(exog_model.predict(full_train[EXOG_FEATURE_COLS]), dtype=float)
    full_test["exog_anom_pred"] = np.asarray(exog_model.predict(full_test[EXOG_FEATURE_COLS]), dtype=float)
    full_train["resid_anom"] = full_train["sales_anom"] - full_train["exog_anom_pred"]
    fallback = float(full_train["resid_anom"].median())
    out_parts: list[pd.DataFrame] = []
    mode_counts: dict[str, int] = {}

    test_series_ids = full_test["series_id"].drop_duplicates().tolist()
    total_series = len(test_series_ids)
    for i, series_id in enumerate(test_series_ids, start=1):
        train_series = full_train[full_train["series_id"] == series_id].sort_values("Date")
        test_series = full_test[full_test["series_id"] == series_id].sort_values("Date").copy()
        horizon = test_series.shape[0]

        if train_series.empty:
            pred = np.full(horizon, fallback, dtype=float)
            mode = "global-fallback"
        else:
            pred, mode = forecast_ets_with_fallback(
                train_series["resid_anom"].to_numpy(dtype=float),
                horizon=horizon,
                seasonal_periods=seasonal_periods,
                fallback_value=fallback,
            )
        mode_counts[mode] = mode_counts.get(mode, 0) + 1
        test_series["pred_resid_anom"] = pred
        test_series["pred_anom"] = test_series["exog_anom_pred"].to_numpy(dtype=float) + test_series[
            "pred_resid_anom"
        ].to_numpy(dtype=float)
        test_series["pred_sales"] = test_series["pred_anom"] + test_series["clim_sales"]
        out_parts.append(
            test_series[["Store", "Dept", "Date", "IsHoliday", "pred_sales", "pred_anom", "series_id"]]
        )

        if i % 400 == 0:
            print(f"  processed {i}/{total_series} series...")

    pred_df = pd.concat(out_parts, ignore_index=True) if out_parts else pd.DataFrame()
    return pred_df, mode_counts


def phase5_default_config() -> dict:
    return {
        "train_path": "train_feat.parquet",
        "test_path": "test_feat.parquet",
        "output_dir": "outputs/baselines/ets",
        "n_folds": 4,
        "val_weeks": 13,
        "max_series": None,
        "seasonal_periods": 52,
        "exog_alpha": 1.0,
    }


def run_phase5(
    config: dict | None = None,
    _ensure_output_dir=ensure_output_dir,
    _load_data=load_data,
    _run_cv_phase5=run_cv,
    _fit_full_phase5=fit_full_and_predict_test,
) -> dict[str, Path]:
    cfg = phase5_default_config()
    if config:
        cfg.update(config)

    warnings.filterwarnings("ignore", category=ConvergenceWarning)
    out_dir = _ensure_output_dir(Path(cfg["output_dir"]))
    train_df, test_df = _load_data(Path(cfg["train_path"]), Path(cfg["test_path"]), max_series=cfg.get("max_series"))

    metrics_df, oof_df, cv_mode_counts = _run_cv_phase5(
        train_df=train_df,
        n_folds=int(cfg["n_folds"]),
        val_weeks=int(cfg["val_weeks"]),
        seasonal_periods=int(cfg["seasonal_periods"]),
        exog_alpha=float(cfg["exog_alpha"]),
    )

    metrics_path = out_dir / "metrics.csv"
    oof_path = out_dir / "oof.parquet"
    metrics_df.to_csv(metrics_path, index=False)
    oof_df.to_parquet(oof_path, index=False)

    test_pred, test_mode_counts = _fit_full_phase5(
        train_df=train_df,
        test_df=test_df,
        seasonal_periods=int(cfg["seasonal_periods"]),
        exog_alpha=float(cfg["exog_alpha"]),
    )
    test_path = out_dir / "test.parquet"
    test_pred.to_parquet(test_path, index=False)

    config_payload = {
        **cfg,
        "model": "ets",
        "target": "sales_anom (via exogenous + ETS residual)",
        "exogenous_features": EXOG_FEATURE_COLS,
        "cv_mode_counts": cv_mode_counts,
        "test_mode_counts": test_mode_counts,
    }
    config_path = out_dir / "config.json"
    config_path.write_text(json.dumps(config_payload, indent=2), encoding="utf-8")

    return {"metrics": metrics_path, "oof": oof_path, "test": test_path, "config": config_path}


def tune_phase5(
    base_config: dict,
    search_space: list[dict],
    search_output_dir: str = "outputs/baselines/ets/search",
    _load_data=load_data,
    _run_cv_phase5=run_cv,
) -> tuple[pd.DataFrame, dict]:
    cfg0 = phase5_default_config()
    cfg0.update(base_config)
    warnings.filterwarnings("ignore", category=ConvergenceWarning)
    train_df, _ = _load_data(Path(cfg0["train_path"]), Path(cfg0["test_path"]), max_series=cfg0.get("max_series"))

    rows: list[dict] = []
    best_cfg: dict | None = None
    best_wmae = float("inf")

    for idx, cand in enumerate(search_space, start=1):
        cfg = {**cfg0, **cand}
        metrics_df, _, _ = _run_cv_phase5(
            train_df=train_df,
            n_folds=int(cfg["n_folds"]),
            val_weeks=int(cfg["val_weeks"]),
            seasonal_periods=int(cfg["seasonal_periods"]),
            exog_alpha=float(cfg["exog_alpha"]),
        )
        mean_row = mean_row_from_metrics(metrics_df)
        wmae_val = float(mean_row["wmae"])
        rows.append(
            {
                "candidate": idx,
                **cand,
                "wmae": wmae_val,
                "mae": float(mean_row.get("mae", np.nan)),
                "rmse": float(mean_row.get("rmse", np.nan)),
            }
        )
        if wmae_val < best_wmae:
            best_wmae = wmae_val
            best_cfg = cfg

    results_df = pd.DataFrame(rows).sort_values("wmae").reset_index(drop=True)
    out_dir = Path(search_output_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    results_df.to_csv(out_dir / "search_results.csv", index=False)
    if best_cfg is None:
        raise RuntimeError("Phase 5 tuning did not evaluate any candidates")
    (out_dir / "best_config.json").write_text(json.dumps(best_cfg, indent=2), encoding="utf-8")
    return results_df, best_cfg

2.8 Phase 6 Core Definitions

import json
from pathlib import Path

from sklearn.ensemble import AdaBoostRegressor
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.tree import DecisionTreeRegressor

def make_model(
    n_estimators: int,
    learning_rate: float,
    max_depth: int,
    loss: str,
):
    tree = DecisionTreeRegressor(max_depth=max_depth, random_state=42)
    # Backward-compatible constructor: newer sklearn uses `estimator`,
    # older versions use `base_estimator`.
    try:
        model = AdaBoostRegressor(
            estimator=tree,
            n_estimators=n_estimators,
            learning_rate=learning_rate,
            loss=loss,
            random_state=42,
        )
    except TypeError:
        model = AdaBoostRegressor(
            base_estimator=tree,
            n_estimators=n_estimators,
            learning_rate=learning_rate,
            loss=loss,
            random_state=42,
        )

    return Pipeline(
        steps=[
            ("imputer", SimpleImputer(strategy="median")),
            ("model", model),
        ]
    )


def phase6_default_config() -> dict:
    return {
        "train_path": "train_feat.parquet",
        "test_path": "test_feat.parquet",
        "output_dir": "outputs/baselines/adaboost",
        "n_folds": 4,
        "val_weeks": 13,
        "max_series": None,
        "n_estimators": 300,
        "learning_rate": 0.03,
        "max_depth": 3,
        "loss": "linear",
    }


def run_phase6(
    config: dict | None = None,
    _ensure_output_dir=ensure_output_dir,
    _load_data=load_data,
    _run_cv=run_cv_tabular_baseline_anomaly,
    _fit_full=fit_full_train_and_predict_test_tabular_anomaly,
    _make_model=make_model,
) -> dict[str, Path]:
    cfg = phase6_default_config()
    if config:
        cfg.update(config)

    out_dir = _ensure_output_dir(Path(cfg["output_dir"]))
    train_df, test_df = _load_data(Path(cfg["train_path"]), Path(cfg["test_path"]), max_series=cfg.get("max_series"))

    metrics_df, oof_df = _run_cv(
        train_df,
        model_builder=lambda: _make_model(
            int(cfg["n_estimators"]),
            float(cfg["learning_rate"]),
            int(cfg["max_depth"]),
            str(cfg["loss"]),
        ),
        n_folds=int(cfg["n_folds"]),
        val_weeks=int(cfg["val_weeks"]),
        lag_orders=LAG_ORDERS,
        feature_cols=ANOM_FEATURE_COLS,
    )

    metrics_path = out_dir / "metrics.csv"
    oof_path = out_dir / "oof.parquet"
    metrics_df.to_csv(metrics_path, index=False)
    oof_df.to_parquet(oof_path, index=False)

    test_pred = _fit_full(
        train_df,
        test_df,
        model_builder=lambda: _make_model(
            int(cfg["n_estimators"]),
            float(cfg["learning_rate"]),
            int(cfg["max_depth"]),
            str(cfg["loss"]),
        ),
        lag_orders=LAG_ORDERS,
        feature_cols=ANOM_FEATURE_COLS,
    )
    test_path = out_dir / "test.parquet"
    test_pred.to_parquet(test_path, index=False)

    config_payload = {
        **cfg,
        "model": "adaboost",
        "target": "sales_anom",
        "features": ANOM_FEATURE_COLS,
        "lags": list(LAG_ORDERS),
    }
    config_path = out_dir / "config.json"
    config_path.write_text(json.dumps(config_payload, indent=2), encoding="utf-8")

    return {"metrics": metrics_path, "oof": oof_path, "test": test_path, "config": config_path}


def tune_phase6(
    base_config: dict,
    search_space: list[dict],
    search_output_dir: str = "outputs/baselines/adaboost/search",
    _load_data=load_data,
    _run_cv=run_cv_tabular_baseline_anomaly,
    _make_model=make_model,
) -> tuple[pd.DataFrame, dict]:
    cfg0 = phase6_default_config()
    cfg0.update(base_config)
    train_df, _ = _load_data(Path(cfg0["train_path"]), Path(cfg0["test_path"]), max_series=cfg0.get("max_series"))

    rows: list[dict] = []
    best_cfg: dict | None = None
    best_wmae = float("inf")

    for idx, cand in enumerate(search_space, start=1):
        cfg = {**cfg0, **cand}
        metrics_df, _ = _run_cv(
            train_df,
            model_builder=lambda: _make_model(
                int(cfg["n_estimators"]),
                float(cfg["learning_rate"]),
                int(cfg["max_depth"]),
                str(cfg["loss"]),
            ),
            n_folds=int(cfg["n_folds"]),
            val_weeks=int(cfg["val_weeks"]),
            lag_orders=LAG_ORDERS,
            feature_cols=ANOM_FEATURE_COLS,
        )
        mean_row = mean_row_from_metrics(metrics_df)
        wmae_val = float(mean_row["wmae"])
        rows.append(
            {
                "candidate": idx,
                **cand,
                "wmae": wmae_val,
                "mae": float(mean_row.get("mae", np.nan)),
                "rmse": float(mean_row.get("rmse", np.nan)),
            }
        )
        if wmae_val < best_wmae:
            best_wmae = wmae_val
            best_cfg = cfg

    results_df = pd.DataFrame(rows).sort_values("wmae").reset_index(drop=True)
    out_dir = Path(search_output_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    results_df.to_csv(out_dir / "search_results.csv", index=False)
    if best_cfg is None:
        raise RuntimeError("Phase 6 tuning did not evaluate any candidates")
    (out_dir / "best_config.json").write_text(json.dumps(best_cfg, indent=2), encoding="utf-8")
    return results_df, best_cfg

2.9 Phase 1 Config and Tuning

PHASE1_CFG = phase1_default_config()
PHASE1_SEARCH = [
    {"kernel": "tricube", "bandwidth": 4, "ridge": 1e-4, "coef_clip": 6.0},
    {"kernel": "tricube", "bandwidth": 6, "ridge": 1e-4, "coef_clip": 6.0},
    {"kernel": "tricube", "bandwidth": 8, "ridge": 5e-4, "coef_clip": 5.0},
    {"kernel": "exponential", "bandwidth": 6, "ridge": 1e-4, "coef_clip": 6.0},
]
PHASE1_RESULTS, PHASE1_BEST_CFG = tune_phase1(PHASE1_CFG, PHASE1_SEARCH)
phase1_artifacts = run_phase1(PHASE1_BEST_CFG)
phase1_artifacts

2.10 Phase 1-lite Config and Run

PHASE1_LITE_CFG = phase1_default_config()
PHASE1_LITE_CFG.update(
    {
        "output_dir": "outputs/phase1_local_linear_lite",
        "feature_mode": "lite",
        "use_interactions": False,
    }
)
phase1_lite_artifacts = run_phase1(PHASE1_LITE_CFG)
phase1_lite_artifacts

2.11 Phase 2 Config and Tuning

PHASE2_CFG = phase2_default_config()
PHASE2_SEARCH = [
    {"fourier_order": 2, "max_eval": 1200, "pred_draws": 20},
    {"fourier_order": 3, "max_eval": 1800, "pred_draws": 40},
    {"fourier_order": 4, "max_eval": 1800, "pred_draws": 40},
]
PHASE2_TUNE_OVERRIDES = {
    "n_folds": 2,
    "val_weeks": 8,
    "max_series": 600,
}
PHASE2_RESULTS, PHASE2_BEST_CFG = tune_phase2(PHASE2_CFG, PHASE2_SEARCH, tune_overrides=PHASE2_TUNE_OVERRIDES)
phase2_artifacts = run_phase2(PHASE2_BEST_CFG)
phase2_artifacts

2.12 Phase 3 Config and Tuning

PHASE3_CFG = phase3_default_config()
PHASE3_SEARCH = [
    {"alpha": 0.01, "l1_ratio": 0.1},
    {"alpha": 0.02, "l1_ratio": 0.2},
    {"alpha": 0.05, "l1_ratio": 0.3},
]
PHASE3_RESULTS, PHASE3_BEST_CFG = tune_phase3(PHASE3_CFG, PHASE3_SEARCH)
phase3_artifacts = run_phase3(PHASE3_BEST_CFG)
phase3_artifacts

2.13 Phase 4 Config and Tuning

PHASE4_CFG = phase4_default_config()
PHASE4_SEARCH = [
    {"n_estimators": 80, "max_depth": 14, "min_samples_leaf": 2, "max_features": "sqrt"},
    {"n_estimators": 120, "max_depth": 18, "min_samples_leaf": 2, "max_features": "sqrt"},
    {"n_estimators": 160, "max_depth": 22, "min_samples_leaf": 1, "max_features": "sqrt"},
]
PHASE4_RESULTS, PHASE4_BEST_CFG = tune_phase4(PHASE4_CFG, PHASE4_SEARCH)
phase4_artifacts = run_phase4(PHASE4_BEST_CFG)
phase4_artifacts

2.14 Phase 5 Config and Tuning

PHASE5_CFG = phase5_default_config()
PHASE5_SEARCH = [
    {"seasonal_periods": 52, "exog_alpha": 0.5},
    {"seasonal_periods": 52, "exog_alpha": 1.0},
    {"seasonal_periods": 53, "exog_alpha": 1.0},
]
PHASE5_RESULTS, PHASE5_BEST_CFG = tune_phase5(PHASE5_CFG, PHASE5_SEARCH)
phase5_artifacts = run_phase5(PHASE5_BEST_CFG)
phase5_artifacts

2.15 Phase 6 Config and Tuning

PHASE6_CFG = phase6_default_config()
PHASE6_SEARCH = [
    {"n_estimators": 200, "learning_rate": 0.03, "max_depth": 2, "loss": "linear"},
    {"n_estimators": 300, "learning_rate": 0.03, "max_depth": 3, "loss": "linear"},
    {"n_estimators": 400, "learning_rate": 0.02, "max_depth": 3, "loss": "square"},
]
PHASE6_RESULTS, PHASE6_BEST_CFG = tune_phase6(PHASE6_CFG, PHASE6_SEARCH)
phase6_artifacts = run_phase6(PHASE6_BEST_CFG)
phase6_artifacts

3 Data and Evaluation Setup

from pathlib import Path
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

PANTONE_GREEN = "#00A651"   # Pantone-style green (approx RGB equivalent)
PANTONE_ORANGE = "#FF6A13"  # Pantone-style orange (approx RGB equivalent)

def find_project_root() -> Path:
    candidates = [Path("."), Path(".."), Path("../..")]
    for c in candidates:
        if (c / "outputs" / "metrics" / "local_linear_cv.csv").exists():
            return c.resolve()
    raise FileNotFoundError("Could not locate project root with outputs/metrics/local_linear_cv.csv")

def split_fold_and_mean(cv_df: pd.DataFrame) -> tuple[pd.DataFrame, pd.Series]:
    fold = cv_df[cv_df["fold"].astype(str) != "mean"].copy()
    mean = cv_df[cv_df["fold"].astype(str) == "mean"].iloc[0]
    return fold, mean

def weekly_actual_vs_pred(oof_df: pd.DataFrame) -> pd.DataFrame:
    agg_map: dict[str, tuple[str, str]] = {
        "actual_sales": ("Weekly_Sales", "sum"),
        "pred_sales": ("pred_sales", "sum"),
    }
    if {"pred_sales_lower", "pred_sales_upper"}.issubset(oof_df.columns):
        agg_map["pred_lower"] = ("pred_sales_lower", "sum")
        agg_map["pred_upper"] = ("pred_sales_upper", "sum")
    out = oof_df.groupby("Date", as_index=False).agg(**agg_map).sort_values("Date")
    return out

def plot_weekly_actual_vs_pred(
    oof_df: pd.DataFrame,
    title: str,
    pred_color: str = PANTONE_ORANGE,
    show_uncertainty: bool = False,
) -> None:
    weekly = weekly_actual_vs_pred(oof_df)
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(weekly["Date"], weekly["actual_sales"], color=PANTONE_GREEN, linewidth=2, label="Actual sales")
    ax.plot(weekly["Date"], weekly["pred_sales"], color=pred_color, linewidth=2, label="Predicted sales")
    if show_uncertainty and {"pred_lower", "pred_upper"}.issubset(weekly.columns):
        ax.fill_between(
            weekly["Date"],
            weekly["pred_lower"],
            weekly["pred_upper"],
            color=pred_color,
            alpha=0.2,
            label="Predictive interval",
        )
    ax.set_title(title)
    ax.set_xlabel("Date")
    ax.set_ylabel("Aggregated Weekly Sales")
    ax.grid(alpha=0.25)
    ax.legend()
    plt.show()

def make_param_df(entries: list[tuple[str, object, str]]) -> pd.DataFrame:
    return pd.DataFrame(entries, columns=["Parameter", "Value", "How it works"])

ROOT = find_project_root()
OUTPUTS = ROOT / "outputs"
BASELINES = OUTPUTS / "baselines"

# Phase 1 (full data)
p1_cv = pd.read_csv(OUTPUTS / "metrics" / "local_linear_cv.csv")
p1_oof = pd.read_parquet(OUTPUTS / "preds" / "local_linear_oof.parquet")
p1_cfg = json.loads((OUTPUTS / "metrics" / "local_linear_config.json").read_text(encoding="utf-8"))
p1_lite_cv = pd.read_csv(OUTPUTS / "phase1_local_linear_lite" / "metrics" / "local_linear_cv.csv")
p1_lite_oof = pd.read_parquet(OUTPUTS / "phase1_local_linear_lite" / "preds" / "local_linear_oof.parquet")
p1_lite_cfg = json.loads(
    (OUTPUTS / "phase1_local_linear_lite" / "metrics" / "local_linear_config.json").read_text(encoding="utf-8")
)
grid_results = pd.read_csv(OUTPUTS / "search" / "local_linear_grid_results.csv")
best_grid_cfg = json.loads((OUTPUTS / "search" / "best_local_linear_config.json").read_text(encoding="utf-8"))

# Phase 2 (full data)
p2_cv = pd.read_csv(OUTPUTS / "phase2_structural_ar_full" / "metrics" / "structural_ar_cv.csv")
p2_oof = pd.read_parquet(OUTPUTS / "phase2_structural_ar_full" / "preds" / "structural_ar_oof.parquet")
p2_cfg = json.loads((OUTPUTS / "phase2_structural_ar_full" / "metrics" / "structural_ar_config.json").read_text(encoding="utf-8"))

# Phase 3, 4, 5, 7 baselines
en_cv = pd.read_csv(BASELINES / "elastic_net" / "metrics.csv")
en_oof = pd.read_parquet(BASELINES / "elastic_net" / "oof.parquet")
en_cfg = json.loads((BASELINES / "elastic_net" / "config.json").read_text(encoding="utf-8"))

rf_cv = pd.read_csv(BASELINES / "random_forest" / "metrics.csv")
rf_oof = pd.read_parquet(BASELINES / "random_forest" / "oof.parquet")
rf_cfg = json.loads((BASELINES / "random_forest" / "config.json").read_text(encoding="utf-8"))

ets_cv = pd.read_csv(BASELINES / "ets" / "metrics.csv")
ets_oof = pd.read_parquet(BASELINES / "ets" / "oof.parquet")
ets_cfg = json.loads((BASELINES / "ets" / "config.json").read_text(encoding="utf-8"))

ad_cv = pd.read_csv(BASELINES / "adaboost" / "metrics.csv")
ad_oof = pd.read_parquet(BASELINES / "adaboost" / "oof.parquet")
ad_cfg = json.loads((BASELINES / "adaboost" / "config.json").read_text(encoding="utf-8"))

for df in [p1_oof, p1_lite_oof, p2_oof, en_oof, rf_oof, ets_oof, ad_oof]:
    df["Date"] = pd.to_datetime(df["Date"])

p1_fold, p1_mean = split_fold_and_mean(p1_cv)
p1_lite_fold, p1_lite_mean = split_fold_and_mean(p1_lite_cv)
p2_fold, p2_mean = split_fold_and_mean(p2_cv)
en_fold, en_mean = split_fold_and_mean(en_cv)
rf_fold, rf_mean = split_fold_and_mean(rf_cv)
ets_fold, ets_mean = split_fold_and_mean(ets_cv)
ad_fold, ad_mean = split_fold_and_mean(ad_cv)

artifact_summary = pd.DataFrame(
    [
        {"Artifact": "Phase 1 OOF rows", "Value": int(p1_oof.shape[0])},
        {"Artifact": "Phase 1-lite OOF rows", "Value": int(p1_lite_oof.shape[0])},
        {"Artifact": "Phase 2 OOF rows", "Value": int(p2_oof.shape[0])},
        {"Artifact": "Elastic Net OOF rows", "Value": int(en_oof.shape[0])},
        {"Artifact": "Random Forest OOF rows", "Value": int(rf_oof.shape[0])},
        {"Artifact": "ETS OOF rows", "Value": int(ets_oof.shape[0])},
        {"Artifact": "AdaBoost OOF rows", "Value": int(ad_oof.shape[0])},
    ]
)
artifact_summary
Artifact Value
0 Phase 1 OOF rows 154386
1 Phase 1-lite OOF rows 154386
2 Phase 2 OOF rows 154386
3 Elastic Net OOF rows 154386
4 Random Forest OOF rows 154386
5 ETS OOF rows 154386
6 AdaBoost OOF rows 154386

Validation metric:

\[ \text{WMAE} = \frac{\sum_i w_i \lvert y_i - \hat{y}_i \rvert}{\sum_i w_i}, \quad w_i = \begin{cases} 5, & \text{if holiday} \\ 1, & \text{otherwise} \end{cases} \]

4 Model 1: Local Linear Anomaly (Phase 1)

Local Linear is a weighted regression that learns a different local coefficient vector for each week-of-year neighborhood and each series. It is simple and interpretable, but sensitive to feature scaling, so standardization is applied before fitting.

4.1 Mathematical Formulation

For series (s) and target seasonal week (w):

\[ \hat{\beta}_{s,w} = \arg\min_{\beta} \sum_{i \in s} K_h\!\left(d(\text{week}_i,w)\right) \left(y_i - \beta_0 - x_i^\top \beta\right)^2 + \lambda \|\beta\|_2^2 \]

and recursive prediction:

\[ \hat{y}^{anom}_{t,s}=\beta_0 + x_{t,s}^\top \beta,\qquad \hat{y}_{t,s}=\hat{y}^{anom}_{t,s}+clim_{t,s} \]

4.2 Plain-Language Meaning of Terms

  • series s: one specific Store + Dept time series.
  • week_i and w: historical week index and target seasonal week where local fitting is centered.
  • K_h(d(.)): a weight that gives more importance to rows close in seasonal calendar (e.g., nearby weeks of year).
  • x_i: the input features for a row (lags, holiday flag, exogenous variables).
  • beta_0, beta: intercept and coefficients learned by local weighted regression.
  • lambda: regularization strength that shrinks coefficients to avoid unstable fits.
  • y^{anom}: anomaly target (real sales minus climatology baseline).
  • clim_{t,s}: baseline seasonal level added back to convert anomaly prediction to sales prediction.

4.3 Workflow

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology per series/store/week]
  B --> C[Create anomalies and lag1/lag2]
  C --> D[Impute missing lag/exogenous values]
  D --> E[Standard-scale Phase 1 features]
  E --> F[Forward-chaining CV]
  F --> G[Fit local weighted ridge by series and seasonal week]
  G --> H[Recursive rollout on validation horizon]
  H --> I[OOF metrics and predictions]
  I --> J[Fit full train and forecast test]

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology per series/store/week]
  B --> C[Create anomalies and lag1/lag2]
  C --> D[Impute missing lag/exogenous values]
  D --> E[Standard-scale Phase 1 features]
  E --> F[Forward-chaining CV]
  F --> G[Fit local weighted ridge by series and seasonal week]
  G --> H[Recursive rollout on validation horizon]
  H --> I[OOF metrics and predictions]
  I --> J[Fit full train and forecast test]

4.4 Tuned Parameters and Effect

make_param_df(
    [
        ("kernel", p1_cfg["kernel"], "Kernel shape over seasonal distance; controls local weighting profile."),
        ("bandwidth", p1_cfg["bandwidth"], "Neighborhood width in week-of-year units; larger means smoother seasonal fit."),
        ("min_samples", p1_cfg["min_samples"], "Minimum active samples to fit a local model; guards against unstable fits."),
        ("ridge", p1_cfg["ridge"], "L2 regularization strength in weighted regression."),
        ("coef_clip", p1_cfg["coef_clip"], "Post-fit coefficient clipping to prevent extreme local slopes."),
        ("anom_clip_scale", p1_cfg["anom_clip_scale"], "Prediction clipping band around anomaly quantiles for recursive stability."),
        ("lags", p1_cfg["lags"], "Autoregressive inputs used in recursive forecasting."),
        ("standard_scaler", p1_cfg.get("standard_scaler", True), "Feature standardization before fitting local regressions."),
    ]
)
Parameter Value How it works
0 kernel tricube Kernel shape over seasonal distance; controls ...
1 bandwidth 6 Neighborhood width in week-of-year units; larg...
2 min_samples 16 Minimum active samples to fit a local model; g...
3 ridge 0.0001 L2 regularization strength in weighted regress...
4 coef_clip 6.0 Post-fit coefficient clipping to prevent extre...
5 anom_clip_scale 2.0 Prediction clipping band around anomaly quanti...
6 lags [1, 2] Autoregressive inputs used in recursive foreca...
7 standard_scaler True Feature standardization before fitting local r...
p1_cfg
{'kernel': 'tricube',
 'bandwidth': 6,
 'min_samples': 16,
 'ridge': 0.0001,
 'coef_clip': 6.0,
 'anom_clip_scale': 2.0,
 'lags': [1, 2],
 'use_interactions': False,
 'interaction_cols': [],
 'feature_cols': ['temp_anom',
  'fuel_anom',
  'MarkDown1',
  'MarkDown2',
  'MarkDown3',
  'MarkDown4',
  'MarkDown5',
  'CPI',
  'Unemployment',
  'sales_anom_lag1',
  'sales_anom_lag2',
  'is_holiday_int'],
 'standard_scaler': True,
 'n_folds': 4,
 'val_weeks': 13,
 'max_series': None}
p1_fold[["fold", "train_start", "train_end", "val_start", "val_end", "wmae", "mae", "rmse"]]
fold train_start train_end val_start val_end wmae mae rmse
0 1 2010-02-05 2011-10-28 2011-11-04 2012-01-27 5989.774614 5916.085624 11465.394577
1 2 2010-02-05 2012-01-27 2012-02-03 2012-04-27 10061.944470 10472.387681 20457.695819
2 3 2010-02-05 2012-04-27 2012-05-04 2012-07-27 11956.737489 11956.737489 21420.299578
3 4 2010-02-05 2012-07-27 2012-08-03 2012-10-26 12294.668865 12479.604773 22343.880598
p1_mean[["wmae", "mae", "rmse"]]
wmae     10075.78136
mae     10206.203892
rmse    18921.817643
Name: 4, dtype: object
plot_weekly_actual_vs_pred(p1_oof, "Model 1: Local Linear | OOF Weekly Sales")

5 Model 1-lite: Local Linear Diagnostic Variant

Phase 1-lite keeps the same Local Linear method but removes the extra exogenous variables (MarkDown1-5, CPI, Unemployment) and uses only temp_anom, fuel_anom plus lags/holiday.

p1_lite_cfg
{'kernel': 'tricube',
 'bandwidth': 6,
 'min_samples': 16,
 'ridge': 0.0001,
 'coef_clip': 6.0,
 'anom_clip_scale': 2.0,
 'lags': [1, 2],
 'feature_mode': 'lite',
 'exogenous_features': ['temp_anom', 'fuel_anom'],
 'use_interactions': False,
 'interaction_cols': [],
 'feature_cols': ['temp_anom',
  'fuel_anom',
  'sales_anom_lag1',
  'sales_anom_lag2',
  'is_holiday_int'],
 'standard_scaler': True,
 'n_folds': 4,
 'val_weeks': 13,
 'max_series': None}
p1_lite_fold[["fold", "train_start", "train_end", "val_start", "val_end", "wmae", "mae", "rmse"]]
fold train_start train_end val_start val_end wmae mae rmse
0 1 2010-02-05 2011-10-28 2011-11-04 2012-01-27 2468.567458 2271.802436 4955.291247
1 2 2010-02-05 2012-01-27 2012-02-03 2012-04-27 2199.682469 2196.184633 4485.720535
2 3 2010-02-05 2012-04-27 2012-05-04 2012-07-27 2028.753667 2028.753667 3948.262807
3 4 2010-02-05 2012-07-27 2012-08-03 2012-10-26 1947.742271 1923.235682 3954.874629
p1_lite_mean[["wmae", "mae", "rmse"]]
wmae    2161.186467
mae     2104.994105
rmse    4336.037305
Name: 4, dtype: object
plot_weekly_actual_vs_pred(p1_lite_oof, "Model 1-lite: Local Linear | OOF Weekly Sales")

6 Model 2: Bayesian Structural AR (Phase 2)

Structural AR combines hierarchical intercepts, autoregressive lags, exogenous drivers, and seasonal Fourier terms in one probabilistic model. It is designed to capture both shared structure and series-specific behavior.

6.1 Mathematical Formulation

In z-score space:

\[ y^{*}_{t,s} \sim \mathcal{N}(\mu_{t,s}, \sigma) \]

\[ \mu_{t,s} = \alpha_s + \beta_{lag}^\top lag_{t,s} + \beta_{exog}^\top exog_{t,s} + \beta_h h_{t,s} + \beta_{tr} trend_t + \beta_f^\top fourier_t \]

with hierarchical intercept:

\[ \alpha_s = \alpha_\mu + \alpha_\sigma z_s,\qquad z_s \sim \mathcal{N}(0,1) \]

6.2 Plain-Language Meaning of Terms

  • y*_{t,s}: normalized anomaly target for week t and series s.
  • mu_{t,s}: model mean prediction before adding random noise.
  • alpha_s: series-specific baseline level; each Store + Dept has its own intercept.
  • beta_lag, lag_{t,s}: coefficients and lag features that capture autoregressive behavior.
  • beta_exog, exog_{t,s}: effects of external drivers (temperature, markdown flags, CPI, unemployment, etc.).
  • beta_h h_{t,s}: holiday contribution.
  • beta_tr trend_t: long-run drift over time.
  • beta_f fourier_t: smooth yearly seasonal pattern via sine/cosine terms.
  • sigma: predictive noise level around the mean.

6.3 Workflow

flowchart LR
  A[Read train parquet] --> B[Compute climatology]
  B --> C[Create anomaly target and lag1/lag2]
  C --> D[Build trend and Fourier seasonal terms]
  D --> E[Fill exogenous NA from train medians]
  E --> F[Standardize lag/exogenous/trend variables]
  F --> G[Fit PyMC structural model via MAP]
  G --> H[Recursive probabilistic rollout by date]
  H --> I[Aggregate draws to mean/sd and intervals]
  I --> J[CV metrics and OOF outputs]

flowchart LR
  A[Read train parquet] --> B[Compute climatology]
  B --> C[Create anomaly target and lag1/lag2]
  C --> D[Build trend and Fourier seasonal terms]
  D --> E[Fill exogenous NA from train medians]
  E --> F[Standardize lag/exogenous/trend variables]
  F --> G[Fit PyMC structural model via MAP]
  G --> H[Recursive probabilistic rollout by date]
  H --> I[Aggregate draws to mean/sd and intervals]
  I --> J[CV metrics and OOF outputs]

6.4 Tuned Parameters and Effect

make_param_df(
    [
        ("max_eval", p2_cfg["max_eval"], "Maximum evaluations in MAP optimization; higher allows deeper convergence."),
        ("pred_draws", p2_cfg["pred_draws"], "Number of stochastic recursive draws for predictive mean/uncertainty."),
        ("fourier_order", p2_cfg["fourier_order"], "Number of sine/cosine seasonal harmonics."),
        ("lag_orders", p2_cfg["lag_orders"], "Autoregressive lags used as structural predictors."),
        ("sigma_clusters", p2_cfg["sigma_clusters"], "Number of heteroskedastic noise clusters; 0 means shared sigma."),
        ("exogenous_features", p2_cfg.get("exogenous_features", []), "External covariates entering the structural mean."),
    ]
)
Parameter Value How it works
0 max_eval 1800 Maximum evaluations in MAP optimization; highe...
1 pred_draws 40 Number of stochastic recursive draws for predi...
2 fourier_order 3 Number of sine/cosine seasonal harmonics.
3 lag_orders [1, 2] Autoregressive lags used as structural predict...
4 sigma_clusters 0 Number of heteroskedastic noise clusters; 0 me...
5 exogenous_features [temp_anom, fuel_anom, MarkDown1, MarkDown2, M... External covariates entering the structural mean.
p2_cfg
{'n_folds': 4,
 'val_weeks': 13,
 'max_eval': 1800,
 'random_seed': 8927,
 'pred_draws': 40,
 'fourier_order': 3,
 'lag_orders': [1, 2],
 'exogenous_features': ['temp_anom',
  'fuel_anom',
  'MarkDown1',
  'MarkDown2',
  'MarkDown3',
  'MarkDown4',
  'MarkDown5',
  'CPI',
  'Unemployment'],
 'sigma_clusters': 0,
 'max_series': 3331,
 'note': 'Structural AR model with hierarchical series intercept, trend, Fourier seasonality, and recursive validation.'}
p2_fold[["fold", "train_start", "train_end", "val_start", "val_end", "wmae", "mae", "rmse", "runtime_sec"]]
fold train_start train_end val_start val_end wmae mae rmse runtime_sec
0 1 2010-02-05 2011-10-28 2011-11-04 2012-01-27 2514.047613 2354.646140 4935.317829 64.718200
1 2 2010-02-05 2012-01-27 2012-02-03 2012-04-27 2192.877850 2207.494598 4400.426532 68.368072
2 3 2010-02-05 2012-04-27 2012-05-04 2012-07-27 2138.504577 2138.504577 3799.486079 68.047948
3 4 2010-02-05 2012-07-27 2012-08-03 2012-10-26 1904.502883 1885.828366 3620.883291 65.612177
p2_mean[["wmae", "mae", "rmse", "runtime_sec"]]
wmae           2187.483231
mae             2146.61842
rmse           4189.028433
runtime_sec      66.686599
Name: 4, dtype: object
plot_weekly_actual_vs_pred(
    p2_oof,
    "Model 2: Structural AR | OOF Weekly Sales with Predictive Uncertainty",
    show_uncertainty=True,
)

7 Model 3: Elastic Net Anomaly Baseline (Phase 3)

Elastic Net is a linear model with combined L1/L2 regularization. It is a robust baseline for correlated tabular features and keeps an interpretable global linear structure.

7.1 Mathematical Formulation

\[ \hat{\beta} = \arg\min_{\beta} \frac{1}{2n}\|y - X\beta\|_2^2 + \alpha\left( \frac{1-l1\_ratio}{2}\|\beta\|_2^2 + l1\_ratio\|\beta\|_1 \right) \]

7.2 Plain-Language Meaning of Terms

  • X beta: linear combination of all features.
  • ||y - X beta||^2: fit error term (how far predictions are from truth).
  • alpha: total regularization strength.
  • l1_ratio: balance between:
  • L1 penalty: pushes less-useful coefficients to exactly zero (feature selection behavior).
  • L2 penalty: smoothly shrinks coefficients for stability under correlated features.

7.3 Workflow

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomalies]
  B --> C[Create lag1/lag2 and calendar features]
  C --> D[Build exogenous feature matrix]
  D --> E[Median imputation + standard scaling]
  E --> F[Fit Elastic Net on train folds]
  F --> G[Recursive forecasting on validation horizon]
  G --> H[OOF metrics and predictions]
  H --> I[Fit full train and forecast test]

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomalies]
  B --> C[Create lag1/lag2 and calendar features]
  C --> D[Build exogenous feature matrix]
  D --> E[Median imputation + standard scaling]
  E --> F[Fit Elastic Net on train folds]
  F --> G[Recursive forecasting on validation horizon]
  G --> H[OOF metrics and predictions]
  H --> I[Fit full train and forecast test]

7.4 Tuned Parameters and Effect

make_param_df(
    [
        ("alpha", en_cfg["alpha"], "Overall regularization strength; larger values shrink coefficients more."),
        ("l1_ratio", en_cfg["l1_ratio"], "Mix between L1 (sparsity) and L2 (stability)."),
        ("max_iter", en_cfg["max_iter"], "Maximum coordinate-descent iterations for convergence."),
        ("lags", en_cfg["lags"], "Lag features used in recursive setup."),
        ("features", en_cfg["features"], "Full tabular interface for model fitting."),
    ]
)
Parameter Value How it works
0 alpha 0.02 Overall regularization strength; larger values...
1 l1_ratio 0.2 Mix between L1 (sparsity) and L2 (stability).
2 max_iter 10000 Maximum coordinate-descent iterations for conv...
3 lags [1, 2] Lag features used in recursive setup.
4 features [lag1, lag2, week_of_year, month, is_holiday_i... Full tabular interface for model fitting.
en_cfg
{'model': 'elastic_net',
 'target': 'sales_anom',
 'features': ['lag1',
  'lag2',
  'week_of_year',
  'month',
  'is_holiday_int',
  'temp_anom',
  'fuel_anom',
  'MarkDown1',
  'MarkDown2',
  'MarkDown3',
  'MarkDown4',
  'MarkDown5',
  'CPI',
  'Unemployment'],
 'lags': [1, 2],
 'alpha': 0.02,
 'l1_ratio': 0.2,
 'max_iter': 10000,
 'n_folds': 4,
 'val_weeks': 13,
 'max_series': None,
 'train_path': 'train_feat.parquet',
 'test_path': 'test_feat.parquet'}
en_fold[["fold", "train_start", "train_end", "val_start", "val_end", "wmae", "mae", "rmse", "runtime_sec"]]
fold train_start train_end val_start val_end wmae mae rmse runtime_sec
0 1 2010-02-05 2011-10-28 2011-11-04 2012-01-27 2472.678568 2292.791482 4919.584363 6.225060
1 2 2010-02-05 2012-01-27 2012-02-03 2012-04-27 2163.772249 2175.479837 4420.761374 8.105888
2 3 2010-02-05 2012-04-27 2012-05-04 2012-07-27 2058.784472 2058.784472 3793.832714 7.717091
3 4 2010-02-05 2012-07-27 2012-08-03 2012-10-26 1833.856224 1803.861985 3645.865425 12.398858
en_mean[["wmae", "mae", "rmse", "runtime_sec"]]
wmae           2132.272878
mae            2082.729444
rmse           4195.010969
runtime_sec       8.611724
Name: 4, dtype: object
plot_weekly_actual_vs_pred(en_oof, "Model 3: Elastic Net | OOF Weekly Sales")

8 Model 4: Random Forest Anomaly Baseline (Phase 4)

Random Forest is an ensemble of decision trees trained on bootstrap samples. It captures nonlinear interactions with little feature engineering and is robust to mixed feature scales.

8.1 Mathematical Formulation

For (B) trees:

\[ \hat{y}(x)=\frac{1}{B}\sum_{b=1}^{B} T_b(x) \]

where each (T_b) is grown on a bootstrap sample and split candidates are randomized via max_features.

8.2 Plain-Language Meaning of Terms

  • T_b(x): prediction from tree b.
  • B: number of trees in the forest.
  • Final prediction is an average of trees, which reduces variance and improves robustness.
  • Bootstrap sampling means each tree sees a slightly different sampled training set.
  • Random split-feature selection (max_features) decorrelates trees, improving ensemble quality.

8.3 Workflow

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomalies]
  B --> C[Create lag1/lag2 + calendar + exogenous features]
  C --> D[Median imputation]
  D --> E[Train Random Forest on fold train]
  E --> F[Recursive multi-step rollout on fold validation]
  F --> G[OOF metrics and predictions]
  G --> H[Refit on full train and forecast test]

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomalies]
  B --> C[Create lag1/lag2 + calendar + exogenous features]
  C --> D[Median imputation]
  D --> E[Train Random Forest on fold train]
  E --> F[Recursive multi-step rollout on fold validation]
  F --> G[OOF metrics and predictions]
  G --> H[Refit on full train and forecast test]

8.4 Tuned Parameters and Effect

make_param_df(
    [
        ("n_estimators", rf_cfg["n_estimators"], "Number of trees; larger reduces variance but increases runtime."),
        ("max_depth", rf_cfg["max_depth"], "Maximum tree depth; controls model complexity."),
        ("min_samples_leaf", rf_cfg["min_samples_leaf"], "Minimum samples per leaf; regularizes tree partitions."),
        ("max_features", rf_cfg["max_features"], "Feature subset size considered at each split."),
        ("lags", rf_cfg["lags"], "Autoregressive lag features used recursively."),
    ]
)
Parameter Value How it works
0 n_estimators 120 Number of trees; larger reduces variance but i...
1 max_depth 18 Maximum tree depth; controls model complexity.
2 min_samples_leaf 2 Minimum samples per leaf; regularizes tree par...
3 max_features sqrt Feature subset size considered at each split.
4 lags [1, 2] Autoregressive lag features used recursively.
rf_cfg
{'model': 'random_forest',
 'target': 'sales_anom',
 'features': ['lag1',
  'lag2',
  'week_of_year',
  'month',
  'is_holiday_int',
  'temp_anom',
  'fuel_anom',
  'MarkDown1',
  'MarkDown2',
  'MarkDown3',
  'MarkDown4',
  'MarkDown5',
  'CPI',
  'Unemployment'],
 'lags': [1, 2],
 'n_estimators': 120,
 'max_depth': 18,
 'min_samples_leaf': 2,
 'max_features': 'sqrt',
 'n_folds': 4,
 'val_weeks': 13,
 'max_series': None,
 'train_path': 'train_feat.parquet',
 'test_path': 'test_feat.parquet'}
rf_fold[["fold", "train_start", "train_end", "val_start", "val_end", "wmae", "mae", "rmse", "runtime_sec"]]
fold train_start train_end val_start val_end wmae mae rmse runtime_sec
0 1 2010-02-05 2011-10-28 2011-11-04 2012-01-27 2443.770022 2243.909149 4898.711983 14.831643
1 2 2010-02-05 2012-01-27 2012-02-03 2012-04-27 2132.907277 2148.134375 4393.405264 20.781198
2 3 2010-02-05 2012-04-27 2012-05-04 2012-07-27 1915.343992 1915.343992 3675.397679 19.938103
3 4 2010-02-05 2012-07-27 2012-08-03 2012-10-26 1738.982368 1719.540629 3554.458799 19.019421
rf_mean[["wmae", "mae", "rmse", "runtime_sec"]]
wmae           2057.750915
mae            2006.732036
rmse           4130.493432
runtime_sec      18.642591
Name: 4, dtype: object
plot_weekly_actual_vs_pred(rf_oof, "Model 4: Random Forest | OOF Weekly Sales")

9 Model 5: ETS Anomaly Baseline (Phase 5)

ETS here is used on residual anomalies after removing an exogenous linear component. This hybrid keeps classical exponential smoothing dynamics while allowing external regressors to explain systematic variation.

9.1 Mathematical Formulation

Exogenous residual decomposition:

\[ y^{anom}_{t,s}=g(x_{t,s}) + r_{t,s},\qquad g(\cdot)\text{ from Ridge regression} \]

Residual component (r_{t,s}) is modeled by ETS:

\[ r_{t,s} = \ell_{t-1,s}+b_{t-1,s}+s_{t-m,s}+\varepsilon_{t,s} \]

and final forecast:

\[ \hat{y}_{t,s}=clim_{t,s}+\hat{g}(x_{t,s})+\hat{r}_{t,s} \]

9.2 Plain-Language Explanation (What ETS Is)

ETS stands for Error, Trend, Seasonality. It is a classical time-series model that keeps three internal states and updates them week by week:

  • level: the current baseline sales level for the series.
  • trend: the direction/slope (increasing or decreasing behavior).
  • seasonality: repeating seasonal pattern (here around yearly weekly cycle).

In this project we first predict the anomaly using exogenous variables (g(x)), then ETS models the remaining residual pattern (r). The final prediction is:

  • exogenous part
  • plus ETS residual dynamics
  • plus climatology baseline to return to sales scale.

9.3 Plain-Language Meaning of Terms

  • y^{anom}: anomaly sales target.
  • g(x): exogenous linear predictor (Ridge) from external variables.
  • r_{t,s}: residual anomaly after removing exogenous part.
  • ell: level state.
  • b: trend state.
  • s: seasonal state.
  • epsilon: random error term.

9.4 Workflow

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomalies]
  B --> C[Fit exogenous Ridge on anomaly target]
  C --> D[Compute residual anomaly = target - exogenous prediction]
  D --> E[Fit per-series ETS on residuals with fallback candidates]
  E --> F[Recursive fold forecasting by series]
  F --> G[Add exogenous component + climatology back]
  G --> H[OOF metrics and test forecasts]

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomalies]
  B --> C[Fit exogenous Ridge on anomaly target]
  C --> D[Compute residual anomaly = target - exogenous prediction]
  D --> E[Fit per-series ETS on residuals with fallback candidates]
  E --> F[Recursive fold forecasting by series]
  F --> G[Add exogenous component + climatology back]
  G --> H[OOF metrics and test forecasts]

9.5 Tuned Parameters and Effect

make_param_df(
    [
        ("seasonal_periods", ets_cfg["seasonal_periods"], "Season length used by ETS seasonal state."),
        ("exog_alpha", ets_cfg.get("exog_alpha", None), "Ridge penalty for exogenous linear component."),
        ("exogenous_features", ets_cfg.get("exogenous_features", []), "Regressors used before ETS residual modeling."),
        ("n_folds", ets_cfg["n_folds"], "Forward-chaining fold count."),
        ("val_weeks", ets_cfg["val_weeks"], "Validation horizon length per fold."),
    ]
)
Parameter Value How it works
0 seasonal_periods 52 Season length used by ETS seasonal state.
1 exog_alpha 1.0 Ridge penalty for exogenous linear component.
2 exogenous_features [temp_anom, fuel_anom, MarkDown1, MarkDown2, M... Regressors used before ETS residual modeling.
3 n_folds 4 Forward-chaining fold count.
4 val_weeks 13 Validation horizon length per fold.
ets_cfg
{'model': 'ets',
 'target': 'sales_anom (via exogenous + ETS residual)',
 'seasonal_periods': 52,
 'exogenous_features': ['temp_anom',
  'fuel_anom',
  'MarkDown1',
  'MarkDown2',
  'MarkDown3',
  'MarkDown4',
  'MarkDown5',
  'CPI',
  'Unemployment'],
 'exog_alpha': 1.0,
 'n_folds': 4,
 'val_weeks': 13,
 'max_series': None,
 'cv_mode_counts': {'level': 10599,
  'level-trend': 1796,
  'global-fallback': 77,
  'last-value-fallback': 51,
  'level-seasonal': 4},
 'test_mode_counts': {'level': 2758,
  'level-trend': 372,
  'level-seasonal': 8,
  'global-fallback': 11,
  'last-value-fallback': 20},
 'train_path': 'train_feat.parquet',
 'test_path': 'test_feat.parquet'}
ets_fold[["fold", "train_start", "train_end", "val_start", "val_end", "wmae", "mae", "rmse", "runtime_sec"]]
fold train_start train_end val_start val_end wmae mae rmse runtime_sec
0 1 2010-02-05 2011-10-28 2011-11-04 2012-01-27 2350.860561 2193.524621 4865.080883 137.155356
1 2 2010-02-05 2012-01-27 2012-02-03 2012-04-27 1965.952051 1969.586155 4035.608919 839.189260
2 3 2010-02-05 2012-04-27 2012-05-04 2012-07-27 1585.766265 1585.766265 3052.148316 952.329147
3 4 2010-02-05 2012-07-27 2012-08-03 2012-10-26 1436.435312 1435.277361 2961.666930 974.137173
ets_mean[["wmae", "mae", "rmse", "runtime_sec"]]
wmae           1834.753547
mae              1796.0386
rmse           3728.626262
runtime_sec     725.702734
Name: 4, dtype: object
plot_weekly_actual_vs_pred(ets_oof, "Model 5: ETS | OOF Weekly Sales")

10 Model 6: AdaBoost Anomaly Baseline (Phase 7)

AdaBoost builds an additive ensemble of shallow trees, where each stage focuses more on difficult residual patterns from previous stages.

10.1 Mathematical Formulation

With base learners (h_m):

\[ F_M(x)=\sum_{m=1}^{M}\nu_m h_m(x) \]

where stage weights depend on chosen boosting loss and learning rate. In practice here, (h_m) are depth-limited regression trees.

10.2 Plain-Language Meaning of Terms

  • h_m(x): base learner at stage m (small regression tree).
  • nu_m: effective contribution of stage m to the final prediction.
  • M: number of boosting stages.
  • learning_rate: shrinkage factor that controls how aggressively each stage updates the model.
  • Boosting logic: each new tree focuses more on errors made by previous trees.

10.3 Workflow

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomaly target]
  B --> C[Create lag1/lag2 + calendar + exogenous features]
  C --> D[Median imputation]
  D --> E[Train AdaBoost regressor on fold train]
  E --> F[Recursive validation rollout]
  F --> G[OOF metrics and predictions]
  G --> H[Refit on full train and produce test forecasts]

flowchart LR
  A[Read train/test parquet] --> B[Compute climatology and anomaly target]
  B --> C[Create lag1/lag2 + calendar + exogenous features]
  C --> D[Median imputation]
  D --> E[Train AdaBoost regressor on fold train]
  E --> F[Recursive validation rollout]
  F --> G[OOF metrics and predictions]
  G --> H[Refit on full train and produce test forecasts]

10.4 Tuned Parameters and Effect

make_param_df(
    [
        ("n_estimators", ad_cfg["n_estimators"], "Number of boosting stages."),
        ("learning_rate", ad_cfg["learning_rate"], "Shrinkage per stage; smaller needs more estimators."),
        ("max_depth", ad_cfg["max_depth"], "Depth of each base regression tree."),
        ("loss", ad_cfg["loss"], "Boosting loss that controls how hard examples are re-weighted."),
        ("lags", ad_cfg["lags"], "Autoregressive lag inputs for recursive forecasting."),
    ]
)
Parameter Value How it works
0 n_estimators 300 Number of boosting stages.
1 learning_rate 0.03 Shrinkage per stage; smaller needs more estima...
2 max_depth 3 Depth of each base regression tree.
3 loss linear Boosting loss that controls how hard examples ...
4 lags [1, 2] Autoregressive lag inputs for recursive foreca...
ad_cfg
{'model': 'adaboost',
 'target': 'sales_anom',
 'features': ['lag1',
  'lag2',
  'week_of_year',
  'month',
  'is_holiday_int',
  'temp_anom',
  'fuel_anom',
  'MarkDown1',
  'MarkDown2',
  'MarkDown3',
  'MarkDown4',
  'MarkDown5',
  'CPI',
  'Unemployment'],
 'lags': [1, 2],
 'n_estimators': 300,
 'learning_rate': 0.03,
 'max_depth': 3,
 'loss': 'linear',
 'n_folds': 4,
 'val_weeks': 13,
 'max_series': None,
 'train_path': 'train_feat.parquet',
 'test_path': 'test_feat.parquet'}
ad_fold[["fold", "train_start", "train_end", "val_start", "val_end", "wmae", "mae", "rmse", "runtime_sec"]]
fold train_start train_end val_start val_end wmae mae rmse runtime_sec
0 1 2010-02-05 2011-10-28 2011-11-04 2012-01-27 2479.961992 2278.508843 4905.691398 144.954249
1 2 2010-02-05 2012-01-27 2012-02-03 2012-04-27 2197.511715 2189.320166 4409.023648 189.685843
2 3 2010-02-05 2012-04-27 2012-05-04 2012-07-27 2014.048705 2014.048705 3866.751795 214.710872
3 4 2010-02-05 2012-07-27 2012-08-03 2012-10-26 1865.985156 1843.294438 3739.207824 239.991155
ad_mean[["wmae", "mae", "rmse", "runtime_sec"]]
wmae           2139.376892
mae            2081.293038
rmse           4230.168666
runtime_sec      197.33553
Name: 4, dtype: object
plot_weekly_actual_vs_pred(ad_oof, "Model 6: AdaBoost | OOF Weekly Sales")

11 Final Comparative

11.1 Full-Data Comparison (Phases 1, 2, 3, 4, 5, 7 + Phase 1-lite Diagnostic)

full_comparison = pd.DataFrame(
    [
        {
            "Model": "Local Linear Anomaly (Phase 1)",
            "Mean WMAE": float(p1_mean["wmae"]),
            "Mean MAE": float(p1_mean["mae"]),
            "Mean RMSE": float(p1_mean["rmse"]),
        },
        {
            "Model": "Structural AR (Phase 2)",
            "Mean WMAE": float(p2_mean["wmae"]),
            "Mean MAE": float(p2_mean["mae"]),
            "Mean RMSE": float(p2_mean["rmse"]),
        },
        {
            "Model": "Local Linear Lite (Phase 1-lite)",
            "Mean WMAE": float(p1_lite_mean["wmae"]),
            "Mean MAE": float(p1_lite_mean["mae"]),
            "Mean RMSE": float(p1_lite_mean["rmse"]),
        },
        {
            "Model": "Elastic Net (Phase 3)",
            "Mean WMAE": float(en_mean["wmae"]),
            "Mean MAE": float(en_mean["mae"]),
            "Mean RMSE": float(en_mean["rmse"]),
        },
        {
            "Model": "Random Forest (Phase 4)",
            "Mean WMAE": float(rf_mean["wmae"]),
            "Mean MAE": float(rf_mean["mae"]),
            "Mean RMSE": float(rf_mean["rmse"]),
        },
        {
            "Model": "ETS (Phase 5)",
            "Mean WMAE": float(ets_mean["wmae"]),
            "Mean MAE": float(ets_mean["mae"]),
            "Mean RMSE": float(ets_mean["rmse"]),
        },
        {
            "Model": "AdaBoost (Phase 7)",
            "Mean WMAE": float(ad_mean["wmae"]),
            "Mean MAE": float(ad_mean["mae"]),
            "Mean RMSE": float(ad_mean["rmse"]),
        },
    ]
).sort_values("Mean WMAE", ascending=True).reset_index(drop=True)
full_comparison
Model Mean WMAE Mean MAE Mean RMSE
0 ETS (Phase 5) 1834.753547 1796.038600 3728.626262
1 Random Forest (Phase 4) 2057.750915 2006.732036 4130.493432
2 Elastic Net (Phase 3) 2132.272878 2082.729444 4195.010969
3 AdaBoost (Phase 7) 2139.376892 2081.293038 4230.168666
4 Local Linear Lite (Phase 1-lite) 2161.186467 2104.994105 4336.037305
5 Structural AR (Phase 2) 2187.483231 2146.618420 4189.028433
6 Local Linear Anomaly (Phase 1) 10075.781360 10206.203892 18921.817643
fig, ax = plt.subplots(figsize=(8, 4))
color_map = {
    "Local Linear Anomaly (Phase 1)": "seagreen",
    "Local Linear Lite (Phase 1-lite)": "darkgreen",
    "Structural AR (Phase 2)": "firebrick",
    "Elastic Net (Phase 3)": "mediumorchid",
    "Random Forest (Phase 4)": "royalblue",
    "ETS (Phase 5)": "teal",
    "AdaBoost (Phase 7)": "darkorange",
}
bar_colors = [color_map[m] for m in full_comparison["Model"]]
ax.bar(
    full_comparison["Model"],
    full_comparison["Mean WMAE"],
    color=bar_colors,
)
ax.set_title("Full-Data Mean WMAE (Lower is Better)")
ax.set_ylabel("Mean WMAE")
ax.grid(axis="y", alpha=0.25)
plt.xticks(rotation=10, ha="right")
plt.tight_layout()
plt.show()

12 Conclusion

  • Full-data comparison is now apples-to-apples for all six models.
  • All tabular baselines share the same anomaly feature interface with lag1 and lag2.
  • The ranking table above is the reference for selecting the best model in this run.